diff --git a/input/methods b/input/methods index c48703f..e7624ab 100644 --- a/input/methods +++ b/input/methods @@ -11,7 +11,7 @@ # RPA RPAx ppRPA F F F # G0F2 evGF2 G0F3 evGF3 - T F F F + F T F F # G0W0 evGW qsGW F F F # G0T0 evGT qsGT diff --git a/src/QuAcK/G0F2.f90 b/src/QuAcK/G0F2.f90 index 3d0905a..0bd5f46 100644 --- a/src/QuAcK/G0F2.f90 +++ b/src/QuAcK/G0F2.f90 @@ -66,6 +66,7 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC, eps = eHF(p) + eHF(a) - eHF(i) - eHF(j) V = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) + Sig(p) = Sig(p) + V*eps/(eps**2 + eta**2) Z(p) = Z(p) - V*(eps**2 - eta**2)/(eps**2 + eta**2)**2 @@ -81,6 +82,7 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC, eps = eHF(p) + eHF(i) - eHF(a) - eHF(b) V = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) + Sig(p) = Sig(p) + V*eps/(eps**2 + eta**2) Z(p) = Z(p) - V*(eps**2 - eta**2)/(eps**2 + eta**2)**2 @@ -111,51 +113,6 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC, call BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE) - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy (singlet) =',EcBSE(1) - write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy (triplet) =',EcBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy =',EcBSE(1) + EcBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F 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,BSE,singlet_manifold,triplet_manifold,eta, & -! nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcAC) - -! if(exchange_kernel) then -! -! EcAC(1) = 0.5d0*EcAC(1) -! EcAC(2) = 1.5d0*EcAC(1) -! -! end if - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! end if - end if + end subroutine G0F2 diff --git a/src/QuAcK/evGF2.f90 b/src/QuAcK/evGF2.f90 index fa99903..c15a205 100644 --- a/src/QuAcK/evGF2.f90 +++ b/src/QuAcK/evGF2.f90 @@ -1,5 +1,5 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold,linearize, & - eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,V,eHF) + eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Perform eigenvalue self-consistent second-order Green function calculation @@ -26,12 +26,14 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: V(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables integer :: nSCF integer :: n_diis + double precision :: EcBSE(nspin) + double precision :: num double precision :: eps double precision :: Conv double precision :: rcond @@ -75,6 +77,7 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol ! Frequency-dependent second-order contribution Sig(:) = 0d0 + Z(:) = 0d0 do p=nC+1,nBas-nR do i=nC+1,nO @@ -82,9 +85,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol do a=nO+1,nBas-nR eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j) + num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) - Sig(p) = Sig(p) & - + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)*eps/(eps**2 + eta**2) + Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do @@ -97,43 +101,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol do b=nO+1,nBas-nR eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b) + num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) - Sig(p) = Sig(p) & - + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)*eps/(eps**2 + eta**2) - - end do - end do - end do - end do - - ! Compute the renormalization factor - - Z(:) = 0d0 - - do p=nC+1,nBas-nR - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nBas-nR - - eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j) - - Z(p) = Z(p) & - - (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)*(eps**2 - eta**2)/(eps**2 + eta**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 = eGF2(p) + eHF(i) - eHF(a) - eHF(b) - - Z(p) = Z(p) & - - (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do @@ -156,7 +127,7 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol ! Print results - call print_evGF2(nBas,nO,nSCF,Conv,eHF,eGF2) + call print_evGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF2) ! DIIS extrapolation @@ -188,4 +159,12 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol end if +! Perform BSE2 calculation + + if(BSE) then + + call BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE) + + end if + end subroutine evGF2