diff --git a/input/methods b/input/methods index 4523991..72b9588 100644 --- a/input/methods +++ b/input/methods @@ -13,7 +13,7 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW - T F F F T F -# G0T0 evGT qsGT - F F F + F F F F F F +# G0T0 evGT qsGT ehG0T0 + F F F T # * unrestricted version available diff --git a/input/options b/input/options index a766382..9ec4b17 100644 --- a/input/options +++ b/input/options @@ -5,14 +5,14 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T T T T + F T F T T # GF: maxSCF thresh DIIS n_diis lin eta renorm reg 256 0.00001 T 5 T 0.0 0 F # GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W reg - 256 0.00001 T 5 T 0.01 F F T F + 256 0.00001 T 5 T 0.0 F F F F # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg 10 0.00001 T 5 T 0.0 F F # ACFDT: AC Kx XBS F T T # BSE: BSE dBSE dTDA evDyn ppBSE BSE2 - F F T F F F + T T T F F F diff --git a/src/GT/ehG0T0.f90 b/src/GT/ehG0T0.f90 new file mode 100644 index 0000000..d1eddf3 --- /dev/null +++ b/src/GT/ehG0T0.f90 @@ -0,0 +1,259 @@ +subroutine ehG0T0(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) + +! Perform ehG0T0 calculation + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: BSE + logical,intent(in) :: BSE2 + logical,intent(in) :: ppBSE + logical,intent(in) :: TDA_T + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: evDyn + 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) :: ENuc + double precision,intent(in) :: ERHF + 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) + double precision,intent(in) :: Vxc(nBas) + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: PHF(nBas,nBas) + +! Local variables + + logical :: print_W = .true. + integer :: ispin + double precision :: EcRPA + double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) + double precision :: EcppBSE(nspin) + double precision :: EcGM + double precision,allocatable :: SigX(:) + double precision,allocatable :: SigC(:) + double precision,allocatable :: Z(:) + 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(:) + +! Output variables + + double precision :: eGT(nBas) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| One-shot eh G0T0 calculation |' + write(*,*)'************************************************' + write(*,*) + +! Initialization + + EcRPA = 0d0 + +! 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 + +! Spin manifold + + ispin = 2 + +! Memory allocation + + allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & + rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS),eGTlin(nBas)) + +!-------------------! +! Compute screening ! +!-------------------! + + call linear_response(ispin,.false.,TDA_T,eta,nBas,nC,nO,nV,nR,nS,1d0, & + eHF,ERI_MO,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + + if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) + +!--------------------------! +! Compute spectral weights ! +!--------------------------! + + call ehGT_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA) + +!------------------------! +! Compute GW self-energy ! +!------------------------! + + call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) + + if(regularize) then + +! call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) +! call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + + 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) + + end if + +!-----------------------------------! +! Solve the quasi-particle equation ! +!-----------------------------------! + + eGTlin(:) = eHF(:) + Z(:)*(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) + + ! Find all the roots of the QP equation if necessary + + ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin) + + end if + +! Compute the RPA correlation energy + + call linear_response(ispin,.false.,TDA_T,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) + +!--------------! +! Dump results ! +!--------------! + + call print_ehG0T0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) + +! Deallocate memory + +! deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin) + +! Plot stuff + +! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,OmRPA,rho_RPA) + +! 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,eHF,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,A3)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! Compute the BSE correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '-------------------------------------------------------------' +! write(*,*) ' Adiabatic connection version of BSE@G0W0 correlation energy ' +! write(*,*) '-------------------------------------------------------------' +! write(*,*) + +! if(doXBS) then + +! write(*,*) '*** scaled screening version (XBS) ***' +! write(*,*) + +! end if + +! call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +! end if + +! if(ppBSE) then + +! call Bethe_Salpeter_pp(TDA_T,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcppBSE) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcppBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcppBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +! if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eGW) +! if(BSE) call ufXBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,OmRPA,rho_RPA) + +! if(BSE) call XBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) + + +end subroutine diff --git a/src/GT/ehGT_excitation_density.f90 b/src/GT/ehGT_excitation_density.f90 new file mode 100644 index 0000000..cc26650 --- /dev/null +++ b/src/GT/ehGT_excitation_density.f90 @@ -0,0 +1,53 @@ +subroutine ehGT_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR) + +! Compute excitation densities + + implicit none + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nR,nS + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: XpY(nS,nS) + double precision,intent(in) :: XmY(nS,nS) + +! Local variables + + integer :: m,jb,p,q,j,b + +! Output variables + + double precision,intent(out) :: rhoL(nBas,nBas,nS) + double precision,intent(out) :: rhoR(nBas,nBas,nS) + + rhoL(:,:,:) = 0d0 + rhoR(:,:,:) = 0d0 + !$OMP PARALLEL & + !$OMP SHARED(nC,nBas,nR,nO,nS,rhoL,rhoR,ERI,XpY,XmY) & + !$OMP PRIVATE(q,p,jb,m) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + do m=1,nS + + rhoL(p,q,m) = rhoL(p,q,m) & + + ERI(p,j,b,q)*0.5d0*(XpY(m,jb) + XmY(m,jb)) & + + ERI(p,b,j,q)*0.5d0*(XpY(m,jb) - XmY(m,jb)) + + rhoR(p,q,m) = rhoR(p,q,m) & + + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*0.5d0*(XpY(m,jb) + XmY(m,jb)) & + + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*0.5d0*(XpY(m,jb) - XmY(m,jb)) + enddo + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + +end subroutine diff --git a/src/GT/ehGT_renormalization_factor.f90 b/src/GT/ehGT_renormalization_factor.f90 new file mode 100644 index 0000000..bf99c9c --- /dev/null +++ b/src/GT/ehGT_renormalization_factor.f90 @@ -0,0 +1,61 @@ +subroutine ehGT_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,Z) + +! Compute renormalization factor for GW + + 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 + integer,intent(in) :: nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rhoL(nBas,nBas,nS) + double precision,intent(in) :: rhoR(nBas,nBas,nS) + +! Local variables + + integer :: p,i,a,m + double precision :: eps + +! Output variables + + double precision,intent(out) :: Z(nBas) + +! Initialize + + Z(:) = 0d0 + +! Occupied part of the correlation self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + do m=1,nS + eps = e(p) - e(i) + Om(m) + Z(p) = Z(p) - rhoL(p,i,m)*rhoR(p,i,m)*(eps/(eps**2 + eta**2))**2 + end do + end do + end do + +! Virtual part of the correlation self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do m=1,nS + eps = e(p) - e(a) - Om(m) + Z(p) = Z(p) - rhoL(p,a,m)*rhoR(p,a,m)*(eps/(eps**2 + eta**2))**2 + end do + end do + end do + +! Compute renormalization factor from derivative of SigC + + Z(:) = 1d0/(1d0 - Z(:)) + +end subroutine diff --git a/src/GT/ehGT_self_energy_diag.f90 b/src/GT/ehGT_self_energy_diag.f90 new file mode 100644 index 0000000..5722dd0 --- /dev/null +++ b/src/GT/ehGT_self_energy_diag.f90 @@ -0,0 +1,74 @@ +subroutine ehGT_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 + + 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 + integer,intent(in) :: nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rhoL(nBas,nBas,nS) + double precision,intent(in) :: rhoR(nBas,nBas,nS) + +! Local variables + + integer :: i,a,p,q,m + double precision :: eps + +! Output variables + + double precision,intent(out) :: SigC(nBas) + double precision,intent(out) :: EcGM + +! Initialize + + SigC(:) = 0d0 + +!----------------------------- +! GW self-energy +!----------------------------- + + ! Occupied part of the correlation self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + do m=1,nS + eps = e(p) - e(i) + Om(m) + SigC(p) = SigC(p) + rhoL(p,i,m)*rhoR(p,i,m)*eps/(eps**2 + eta**2) + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do m=1,nS + eps = e(p) - e(a) - Om(m) + SigC(p) = SigC(p) + rhoL(p,a,m)*rhoR(p,a,m)*eps/(eps**2 + eta**2) + end do + end do + end do + + ! GM correlation energy + + EcGM = 0d0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do m=1,nS + eps = e(a) - e(i) + Om(m) + EcGM = EcGM - 2d0*rhoL(a,i,m)*rhoR(a,i,m)*eps/(eps**2 + eta**2) + end do + end do + end do + +end subroutine diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 0b82089..4627a3d 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -191,7 +191,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD ! Deallocate memory - deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin) +! deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin) ! Plot stuff @@ -280,4 +280,10 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD end if +! if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eGW) +! if(BSE) call ufXBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,OmRPA,rho_RPA) + + if(BSE) call XBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) + + end subroutine G0W0 diff --git a/src/GW/XBSE.f90 b/src/GW/XBSE.f90 new file mode 100644 index 0000000..9e8eb92 --- /dev/null +++ b/src/GW/XBSE.f90 @@ -0,0 +1,133 @@ +subroutine XBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) + +! Compute the Bethe-Salpeter excitation energies + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: singlet + logical,intent(in) :: triplet + + 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 + integer,intent(in) :: nS + double precision,intent(in) :: eW(nBas) + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: ispin + integer :: isp_W + + double precision :: EcRPA + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) + + double precision,allocatable :: OmBSE(:,:) + double precision,allocatable :: XpY_BSE(:,:,:) + double precision,allocatable :: XmY_BSE(:,:,:) + + double precision,allocatable :: KA_sta(:,:) + double precision,allocatable :: KB_sta(:,:) + + double precision,allocatable :: W(:,:,:,:) + +! Output variables + + double precision,intent(out) :: EcBSE(nspin) + +! Memory allocation + + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & + KA_sta(nS,nS),KB_sta(nS,nS),OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) + +!--------------------------------- +! Compute (singlet) RPA screening +!--------------------------------- + + isp_W = 1 + EcRPA = 0d0 + + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + + call XBSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta,eW,eGW) + call XBSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta,eW) + +!------------------- +! Singlet manifold +!------------------- + + if(singlet) then + + ispin = 1 + EcBSE(ispin) = 0d0 + + ! Compute BSE excitation energies + + call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,KA_sta,KB_sta, & + EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) + call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + + !------------------------------------------------- + ! Compute the dynamical screening at the BSE level + !------------------------------------------------- + + if(dBSE) then + + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eW,dipole_int,OmRPA,rho_RPA, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + end if + + end if + +!------------------- +! Triplet manifold +!------------------- + + if(triplet) then + + ispin = 2 + EcBSE(ispin) = 0d0 + + ! Compute BSE excitation energies + + call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,KA_sta,KB_sta, & + EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + + !------------------------------------------------- + ! Compute the dynamical screening at the BSE level + !------------------------------------------------- + + if(dBSE) then + + ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) + + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eW,dipole_int,OmRPA,rho_RPA, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + + end if + + end if + +end subroutine diff --git a/src/GW/XBSE_static_kernel_KA.f90 b/src/GW/XBSE_static_kernel_KA.f90 new file mode 100644 index 0000000..46167f8 --- /dev/null +++ b/src/GW/XBSE_static_kernel_KA.f90 @@ -0,0 +1,119 @@ +subroutine XBSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,WA,eW,eGW) + +! Compute the BSE static kernel for the resonant block + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: eW(nBas) + double precision,intent(in) :: eGW(nBas) + +! Local variables + + double precision :: chi + double precision :: eps + double precision :: num,den,ei,ea + integer :: i,j,k,a,b,c,ia,jb,m + double precision,external :: Kronecker_delta + +! Output variables + + double precision,intent(out) :: WA(nS,nS) + +! Initialize + + WA(:,:) = 0d0 + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + ! virtual quasiparticle term + + ea = 0d0 + do m=1,nS + do k=nC+1,nO + num = 1d0*rho(a,k,m)*rho(b,k,m) + den = eW(a) - eW(k) + Om(m) + + ea = ea + Kronecker_delta(i,j)*num*den/(den**2+eta**2) + + den = eW(b) - eW(k) + Om(m) + + ea = ea + Kronecker_delta(i,j)*num*den/(den**2+eta**2) + + end do + do c=nO+1,nBas-nR + num = 1d0*rho(a,c,m)*rho(b,c,m) + den = eW(a) - eW(c) - Om(m) + + ea = ea + Kronecker_delta(i,j)*num*den/(den**2+eta**2) + + den = eW(b) - eW(c) - Om(m) + + ea = ea + Kronecker_delta(i,j)*num*den/(den**2+eta**2) + + end do + end do + + ! occupied quasiparticle term + + ei = 0d0 + do m=1,nS + do k=nC+1,nO + num = 1d0*rho(i,k,m)*rho(j,k,m) + den = eW(i) - eW(k) + Om(m) + + ei = ei + Kronecker_delta(a,b)*num*den/(den**2+eta**2) + + den = eW(j) - eW(k) + Om(m) + + ei = ei + Kronecker_delta(a,b)*num*den/(den**2+eta**2) + + end do + do c=nO+1,nBas-nR + num = 1d0*rho(i,c,m)*rho(j,c,m) + den = eW(i) - eW(c) - Om(m) + + ei = ei + Kronecker_delta(a,b)*num*den/(den**2+eta**2) + + den = eW(j) - eW(c) - Om(m) + + ei = ei + Kronecker_delta(a,b)*num*den/(den**2+eta**2) + + end do + end do + + ! kernel term + + chi = 0d0 + do m=1,nS + eps = Om(m)**2 + eta**2 + chi = chi + rho(i,j,m)*rho(a,b,m)*Om(m)/eps!& +! - rho(i,b,m)*rho(a,j,m)*Om(m)/eps + enddo + +! WA(ia,jb) = WA(ia,jb) + lambda*ERI(i,b,j,a) - 4d0*lambda*chi & +! - (eGW(a) - eW(a))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & +! + (eGW(i) - eW(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) + WA(ia,jb) = WA(ia,jb) - ea + ei + lambda*ERI(i,b,j,a) - 4d0*lambda*chi + + enddo + enddo + enddo + enddo + +end subroutine diff --git a/src/GW/XBSE_static_kernel_KB.f90 b/src/GW/XBSE_static_kernel_KB.f90 new file mode 100644 index 0000000..2ca9b65 --- /dev/null +++ b/src/GW/XBSE_static_kernel_KB.f90 @@ -0,0 +1,59 @@ +subroutine XBSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,KB) + +! Compute the BSE static kernel for the coupling block + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: KB(nS,nS) + +! Initialize + + KB(:,:) = 0d0 + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + chi = 0d0 + do kc=1,nS + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps & + - rho(i,a,kc)*rho(j,b,kc)*Omega(kc)/eps + enddo + + KB(ia,jb) = KB(ia,jb) + lambda*ERI(i,j,b,a) - 4d0*lambda*chi + + enddo + enddo + enddo + enddo + +end subroutine diff --git a/src/GW/ufBSE.f90 b/src/GW/ufBSE.f90 index 426b4bc..eaf2fa8 100644 --- a/src/GW/ufBSE.f90 +++ b/src/GW/ufBSE.f90 @@ -25,6 +25,7 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) integer :: i,j,k,l integer :: a,b,c,d integer :: ia,jb,kc,iajb,kcld + integer,parameter :: maxH = 20 double precision :: tmp integer :: n1h1p,n2h2p,nH @@ -66,7 +67,7 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) H(:,:) = 0d0 !---------------------------! -! Compute GW supermatrix ! +! Compute BSE supermatrix ! !---------------------------! ! ! ! | A -Ve -Vh | ! @@ -202,7 +203,7 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) '|','#','|','Omega (eV)','|','Z','|' write(*,*)'-------------------------------------------' - do s=1,nH + do s=1,min(nH,maxH) if(Z(s) > 1d-7) & write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',s,'|',Om(s)*HaToeV,'|',Z(s),'|' diff --git a/src/GW/ufXBSE.f90 b/src/GW/ufXBSE.f90 new file mode 100644 index 0000000..31acdbc --- /dev/null +++ b/src/GW/ufXBSE.f90 @@ -0,0 +1,268 @@ +subroutine ufXBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,OmRPA,sERI) + +! Unfolded BSE+ equations + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: sERI(nBas,nBas,nS) + +! Local variables + + integer :: s + integer :: i,j,k,l + integer :: a,b,c,d + integer :: ia,jb,kc,iajb,kcld + integer,parameter :: maxH = 20 + double precision :: eps1,eps2 + double precision :: Ve,Vh,C2h2p + + integer :: n1h1p,n2h2p,nH + double precision,external :: Kronecker_delta + + double precision,allocatable :: H(:,:) + double precision,allocatable :: X(:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: Z(:) + +! Output variables + +! Hello world + + write(*,*) + write(*,*)'**********************************************' + write(*,*)'| Unfolded BSE+ calculation |' + write(*,*)'**********************************************' + write(*,*) + +! TDA for W + + write(*,*) 'Tamm-Dancoff approximation by default!' + write(*,*) + +! Dimension of the supermatrix + + n1h1p = nO*nV + n2h2p = nO*nO*nV*nV + nH = n1h1p + n2h2p + +! Memory allocation + + allocate(H(nH,nH),X(nH,nH),Om(nH),Z(nH)) + +! Initialization + + H(:,:) = 0d0 + +!---------------------------! +! Compute BSE+ supermatrix ! +!---------------------------! +! ! +! | A Ve-Vh | ! +! H = | | ! +! | Ve-Vh C2h2p | ! +! ! +!---------------------------! + + !---------! + ! Block A ! + !---------! + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + H(ia,jb) = (eHF(a) - eHF(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + + 2d0*ERI(i,b,a,j) - ERI(i,b,j,a) + + do kc=1,nS + + do l=nC+1,nO + eps1 = 1d0/(eHF(a) - eHF(l) + OmRPA(kc)) + eps2 = 1d0/(eHF(b) - eHF(l) + OmRPA(kc)) + H(ia,jb) = H(ia,jb) + Kronecker_delta(i,j)*sERI(a,l,kc)*sERI(b,l,kc)*(eps1+eps2) + enddo + + do d=nO+1,nBas-nR + eps1 = 1d0/(- eHF(i) + eHF(d) + OmRPA(kc)) + eps2 = 1d0/(- eHF(j) + eHF(d) + OmRPA(kc)) + H(ia,jb) = H(ia,jb) + Kronecker_delta(a,b)*sERI(i,d,kc)*sERI(j,d,kc)*(eps1+eps2) + enddo + + eps1 = 1d0/(eHF(a) - eHF(i) + OmRPA(kc)) + eps2 = 1d0/(eHF(b) - eHF(j) + OmRPA(kc)) + H(ia,jb) = H(ia,jb) - 2d0*sERI(i,a,kc)*sERI(j,b,kc)*(eps1+eps2) + + end do + + end do + end do + + end do + end do + + !----------------! + ! Blocks Vp & Ve ! + !----------------! + + iajb=0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do j=nC+1,nO + do b=nO+1,nBas-nR + iajb = iajb + 1 + + kc = 0 + do k=nC+1,nO + do c=nO+1,nBas-nR + kc = kc + 1 + + Ve = sqrt(2d0)*Kronecker_delta(k,j)*ERI(b,a,c,i) + Vh = sqrt(2d0)*Kronecker_delta(b,c)*ERI(a,k,i,j) + + H(n1h1p+iajb,kc ) = Ve - Vh + H(kc ,n1h1p+iajb) = Ve - Vh + + end do + end do + + end do + end do + end do + end do + +! iajb=0 +! ia = 0 +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! ia = ia + 1 +! do j=nC+1,nO +! do b=nO+1,nBas-nR +! iajb = iajb + 1 + +! kc = 0 +! do k=nC+1,nO +! do c=nO+1,nBas-nR +! kc = kc + 1 + +! Ve = sqrt(2d0)*Kronecker_delta(k,j)*sERI(b,c,ia) +! Vh = sqrt(2d0)*Kronecker_delta(b,c)*sERI(k,j,ia) + +! H(n1h1p+iajb,kc ) = Ve - Vh +! H(kc ,n1h1p+iajb) = Ve - Vh +! +! end do +! end do + +! end do +! end do +! end do +! end do + + !-------------! + ! Block 2h2p ! + !-------------! + + iajb = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do j=nC+1,nO + do b=nO+1,nBas-nR + iajb = iajb + 1 + + kcld = 0 + do k=nC+1,nO + do c=nO+1,nBas-nR + do l=nC+1,nO + do d=nO+1,nBas-nR + kcld = kcld + 1 + + C2h2p = ((eHF(a) + eHF(b) - eHF(i) - eHF(j))*Kronecker_delta(i,k)*Kronecker_delta(a,c) & + + 2d0*ERI(a,k,i,c))*Kronecker_delta(j,l)*Kronecker_delta(b,d) + + H(n1h1p+iajb,n1h1p+kcld) = C2h2p + + end do + end do + end do + end do + + end do + end do + end do + end do + +! iajb = 0 +! ia = 0 +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! ia = ia + 1 +! do j=nC+1,nO +! do b=nO+1,nBas-nR +! iajb = iajb + 1 + +! H(n1h1p+iajb,n1h1p+iajb) = Om(ia) + eHF(b) - eHF(j) + +! end do +! end do +! end do +! end do + +!-------------------------! +! Diagonalize supermatrix ! +!-------------------------! + + X(:,:) = H(:,:) + call diagonalize_matrix(nH,X,Om) + +!-----------------! +! Compute weights ! +!-----------------! + + Z(:) = 0d0 + do s=1,nH + do ia=1,n1h1p + Z(s) = Z(s) + X(ia,s)**2 + end do + end do + +!--------------! +! Dump results ! +!--------------! + + write(*,*)'-------------------------------------------' + write(*,*)' BSE+ excitation energies (eV) ' + write(*,*)'-------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') & + '|','#','|','Omega (eV)','|','Z','|' + write(*,*)'-------------------------------------------' + + do s=1,min(nH,maxH) + if(Z(s) > 1d-7) & + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',s,'|',Om(s)*HaToeV,'|',Z(s),'|' + enddo + + write(*,*)'-------------------------------------------' + write(*,*) + +end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 756d59c..b89450c 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -17,6 +17,7 @@ program QuAcK logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW logical :: doG0T0,doevGT,doqsGT + logical :: doehG0T0 integer :: nNuc,nBas,nBasCABS integer :: nEl(nspin) @@ -169,7 +170,8 @@ program QuAcK doG0F3,doevGF3, & doG0W0,doevGW,doqsGW,doSRGqsGW, & doufG0W0,doufGW, & - doG0T0,doevGT,doqsGT) + doG0T0,doevGT,doqsGT, & + doehG0T0) ! Read options for methods @@ -1099,6 +1101,8 @@ program QuAcK write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_ufGW,' seconds' write(*,*) + if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) + end if !------------------------------------------------------------------------ @@ -1109,7 +1113,7 @@ program QuAcK call cpu_time(start_ufGW) call ufGW(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) - call CCGW(maxSCF_CC,thresh_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) +! call CCGW(maxSCF_CC,thresh_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_ufGW) t_ufGW = end_ufGW - start_ufGW @@ -1219,6 +1223,35 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Perform ehG0T0 calculatiom +!------------------------------------------------------------------------ + + eG0T0(:,:) = eHF(:,:) + + if(doehG0T0) then + + call cpu_time(start_G0T0) + + if(unrestricted) then + + print*,'!!! ehG0T0 NYI at the unrestricted level !!!' + + else + + call ehG0T0(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 + + call cpu_time(end_G0T0) + + t_G0T0 = end_G0T0 - start_G0T0 + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Compute FCI !------------------------------------------------------------------------ diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 990cb3c..459516a 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -8,7 +8,8 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doG0F3,doevGF3, & doG0W0,doevGW,doqsGW,doSRGqsGW, & doufG0W0,doufGW, & - doG0T0,doevGT,doqsGT) + doG0T0,doevGT,doqsGT, & + doehG0T0) ! Read desired methods @@ -25,6 +26,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & 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 ! Local variables @@ -159,10 +161,11 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Read GT methods read(1,*) - read(1,*) answer1,answer2,answer3 - if(answer1 == 'T') doG0T0 = .true. - if(answer2 == 'T') doevGT = .true. - if(answer3 == 'T') doqsGT = .true. + 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. ! Close file with geometry specification