diff --git a/input/methods b/input/methods index 4d84f78..5f36a66 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 F F + F T F F F F # G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh F F F F F F # * unrestricted version available diff --git a/input/options b/input/options index 4423bd6..e73a4a6 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability - 999 0.0000001 T 5 1 1 T 0.0 T + 5000 0.0000001 T 5 1 1 T 0.0 T # MP: reg F # CC: maxSCF thresh DIIS n_diis @@ -9,7 +9,7 @@ # 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 TDA_W reg - 256 0.00001 T 5 T 0.00367493 F F + 256 0.00001 T 5 F 0.0 F F # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg 256 0.00001 T 5 F 0.0 F F # ACFDT: AC Kx XBS diff --git a/src/GW/GW_SigC.f90 b/src/GW/GW_SigC.f90 index 78564f4..bff10fc 100644 --- a/src/GW/GW_SigC.f90 +++ b/src/GW/GW_SigC.f90 @@ -32,21 +32,21 @@ double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Occupied part of the correlation self-energy do i=nC+1,nO - do m=1,nS - eps = w - e(i) + Om(m) - num = 2d0*rho(p,i,m)**2 - GW_SigC = GW_SigC + num*eps/(eps**2 + eta**2) - enddo - enddo + do m=1,nS + eps = w - e(i) + Om(m) + num = 2d0*rho(p,i,m)**2 + GW_SigC = GW_SigC + num*eps/(eps**2 + eta**2) + end do + end do ! Virtual part of the correlation self-energy do a=nO+1,nBas-nR - do m=1,nS - eps = w - e(a) - Om(m) - num = 2d0*rho(p,a,m)**2 - GW_SigC = GW_SigC + num*eps/(eps**2 + eta**2) - enddo - enddo + do m=1,nS + eps = w - e(a) - Om(m) + num = 2d0*rho(p,a,m)**2 + GW_SigC = GW_SigC + num*eps/(eps**2 + eta**2) + end do + end do end function diff --git a/src/GW/GW_dSigC.f90 b/src/GW/GW_dSigC.f90 index 19b14c7..aad9921 100644 --- a/src/GW/GW_dSigC.f90 +++ b/src/GW/GW_dSigC.f90 @@ -32,21 +32,21 @@ double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Occupied part of the correlation self-energy do i=nC+1,nO - do m=1,nS - eps = w - e(i) + Om(m) - num = 2d0*rho(p,i,m)**2 - GW_dSigC = GW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo + do m=1,nS + eps = w - e(i) + Om(m) + num = 2d0*rho(p,i,m)**2 + GW_dSigC = GW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + end do ! Virtual part of the correlation self-energy do a=nO+1,nBas-nR - do m=1,nS - eps = w - e(a) - Om(m) - num = 2d0*rho(p,a,m)**2 - GW_dSigC = GW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo + do m=1,nS + eps = w - e(a) - Om(m) + num = 2d0*rho(p,a,m)**2 + GW_dSigC = GW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + end do end function diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 96d234a..75463f6 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -147,9 +147,17 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons ! Find graphical solution of the QP equation + write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + write(*,*) + do is=1,nspin + + write(*,*)'-----------------------------------------------------' + if(is==1) write(*,*)' Spin-up orbitals ' + if(is==2) write(*,*)' Spin-down orbitals ' + call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eHF(:,is), & - Om,rho(:,:,:,is),eHF(:,is),eGW(:,is),Z(:,is)) + Om,rho(:,:,:,is),eGWlin(:,is),eHF(:,is),eGW(:,is),Z(:,is)) end do end if diff --git a/src/GW/UGW_QP_graph.f90 b/src/GW/UGW_QP_graph.f90 index 5c5d9ca..dfde236 100644 --- a/src/GW/UGW_QP_graph.f90 +++ b/src/GW/UGW_QP_graph.f90 @@ -1,4 +1,4 @@ -subroutine UGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z) +subroutine UGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) ! Compute the graphical solution of the QP equation @@ -13,12 +13,14 @@ subroutine UGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z) integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS + double precision,intent(in) :: eta double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS,nspin) double precision,intent(in) :: eGWlin(nBas) + double precision,intent(in) :: eOld(nBas) ! Local variables @@ -27,7 +29,7 @@ subroutine UGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z) integer,parameter :: maxIt = 64 double precision,parameter :: thresh = 1d-6 double precision,external :: UGW_SigC,UGW_dSigC - double precision :: sigC,dsigC + double precision :: SigC,dSigC double precision :: f,df double precision :: w @@ -38,44 +40,40 @@ subroutine UGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z) ! Run Newton's algorithm to find the root - do p=nC+1,nBas-nR + write(*,*)'-----------------------------------------------------' + write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','e_GWlin (eV)','e_GW (eV)','Z' + write(*,*)'-----------------------------------------------------' - write(*,*) '-----------------' - write(*,'(A10,I3)') 'Orbital ',p - write(*,*) '-----------------' + do p=nC+1,nBas-nR w = eGWlin(p) nIt = 0 f = 1d0 - write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f do while (abs(f) > thresh .and. nIt < maxIt) nIt = nIt + 1 - sigC = UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho) - dsigC = UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho) - f = w - eHF(p) - sigC - df = 1d0/(1d0 - dsigC) + SigC = UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + dSigC = UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + f = w - eHF(p) - SigC + df = 1d0/(1d0 - dSigC) w = w - df*f - - write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,df,f end do if(nIt == maxIt) then - write(*,*) 'Newton root search has not converged!' eGW(p) = eGWlin(p) + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,eGWlin(p)*HaToeV,eGW(p)*HaToeV,Z(p),'Cvg Failed!' else eGW(p) = w Z(p) = df - write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(p)*HaToeV - write(*,*) + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,eGWlin(p)*HaToeV,eGW(p)*HaToeV,Z(p) end if diff --git a/src/GW/UGW_SigC.f90 b/src/GW/UGW_SigC.f90 index 1dc5d2d..488c983 100644 --- a/src/GW/UGW_SigC.f90 +++ b/src/GW/UGW_SigC.f90 @@ -22,7 +22,7 @@ double precision function UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Local variables - integer :: i,a,jb + integer :: i,a,m double precision :: num,eps ! Initialize @@ -32,9 +32,9 @@ double precision function UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Occupied part of the correlation self-energy do i=nC+1,nO - do jb=1,nS - eps = w - e(i) + Om(jb) - num = rho(p,i,jb)**2 + do m=1,nS + eps = w - e(i) + Om(m) + num = rho(p,i,m)**2 UGW_SigC = UGW_SigC + num*eps/(eps**2 + eta**2) end do end do @@ -42,9 +42,9 @@ double precision function UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Virtual part of the correlation self-energy do a=nO+1,nBas-nR - do jb=1,nS - eps = w - e(a) - Om(jb) - num = rho(p,a,jb)**2 + do m=1,nS + eps = w - e(a) - Om(m) + num = rho(p,a,m)**2 UGW_SigC = UGW_SigC + num*eps/(eps**2 + eta**2) end do end do diff --git a/src/GW/UGW_dSigC.f90 b/src/GW/UGW_dSigC.f90 index 7e9ac65..8a82fc7 100644 --- a/src/GW/UGW_dSigC.f90 +++ b/src/GW/UGW_dSigC.f90 @@ -22,7 +22,7 @@ double precision function UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Local variables - integer :: i,a,jb + integer :: i,a,m double precision :: num,eps ! Initialize @@ -32,20 +32,20 @@ double precision function UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Occupied part of the correlation self-energy do i=nC+1,nO - do jb=1,nS - eps = w - e(i) + Om(jb) - num = rho(p,i,jb)**2 - UGW_dSigC = UGW_dSigC + num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + do m=1,nS + eps = w - e(i) + Om(m) + num = rho(p,i,m)**2 + UGW_dSigC = UGW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do ! Virtual part of the correlation self-energy do a=nO+1,nBas-nR - do jb=1,nS - eps = w - e(a) - Om(jb) - num = rho(p,a,jb)**2 - UGW_dSigC = UGW_dSigC + num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + do m=1,nS + eps = w - e(a) - Om(m) + num = rho(p,a,m)**2 + UGW_dSigC = UGW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index bcd8c07..e5ab829 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -166,8 +166,13 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, write(*,*) do is=1,nspin + + write(*,*)'-----------------------------------------------------' + if(is==1) write(*,*)' Spin-up orbitals ' + if(is==2) write(*,*)' Spin-down orbitals ' + call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eHF(:,is), & - Om,rho(:,:,:,is),eOld(:,is),eGW(:,is),Z(:,is)) + Om,rho(:,:,:,is),eOld(:,is),eOld(:,is),eGW(:,is),Z(:,is)) end do end if