From 172ca40f370eebe1815f384bae61b689584116f0 Mon Sep 17 00:00:00 2001 From: pfloos Date: Mon, 30 Sep 2024 23:06:00 +0200 Subject: [PATCH] saving worj in ccGW --- src/GW/RGW.f90 | 8 +- src/GW/ccRG0W0.f90 | 146 +++++++++++----------- src/GW/ccRG0W0_TDA.f90 | 263 ++++++++++++++++++++++++++++++++++++++ src/GW/ccRG0W0_mat.f90 | 278 ----------------------------------------- 4 files changed, 337 insertions(+), 358 deletions(-) create mode 100644 src/GW/ccRG0W0_TDA.f90 delete mode 100644 src/GW/ccRG0W0_mat.f90 diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index 0e2a63a..73e723a 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -162,12 +162,13 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii ! Perform CC-based G0W0 calculation !------------------------------------------------------------------------ - doccG0W0 = .false. + doccG0W0 = .true. if(doccG0W0) then call wall_time(start_GW) - call ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,nS,ERI_MO,ENuc,ERHF,eHF) +! call ccRG0W0_TDA(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -186,8 +187,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii if(doccGW) then call wall_time(start_GW) -! call ccRGW(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) -! call ccRGW_mat(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call ccRGW(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/GW/ccRG0W0.f90 b/src/GW/ccRG0W0.f90 index a8b8d54..d3ec39e 100644 --- a/src/GW/ccRG0W0.f90 +++ b/src/GW/ccRG0W0.f90 @@ -1,4 +1,4 @@ -subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) +subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,nS,ERI,ENuc,ERHF,eHF) ! CC-based GW module @@ -17,6 +17,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH 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) :: eHF(nOrb) @@ -27,17 +28,29 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH integer :: p,q,r,s integer :: i,j,k,l integer :: a,b,c,d + integer :: m + integer :: isp_W + logical :: TDA_W + logical :: dRPA integer :: nSCF double precision :: Conv + double precision :: EcRPA double precision,allocatable :: Sig(:) double precision,allocatable :: Z(:) - double precision,allocatable :: del(:,:,:) - double precision,allocatable :: vec(:,:,:) - double precision,allocatable :: res(:,:,:) - double precision,allocatable :: amp(:,:,:) + double precision,allocatable :: del(:,:) + double precision,allocatable :: vec(:,:) + double precision,allocatable :: res(:,:) + double precision,allocatable :: amp(:,:) + + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: rho(:,:,:) integer :: n_diis double precision :: rcond @@ -55,16 +68,39 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH ! Memory allocation - allocate(del(nO,nV,nOrb)) - allocate(vec(nO,nV,nOrb)) - allocate(res(nO,nV,nOrb)) - allocate(amp(nO,nV,nOrb)) + allocate(del(nS,nOrb)) + allocate(vec(nS,nOrb)) + allocate(res(nS,nOrb)) + allocate(amp(nS,nOrb)) allocate(Sig(nOrb)) allocate(Z(nOrb)) - allocate(r_diis(nO*nV*nOrb,max_diis)) - allocate(t_diis(nO*nV*nOrb,max_diis)) + allocate(r_diis(nS*nOrb,max_diis)) + allocate(t_diis(nS*nOrb,max_diis)) + +!-------------------! +! Compute screening ! +!-------------------! + + ! Spin manifold + + isp_W = 1 + TDA_W = .false. + dRPA = .true. + + ! Memory allocation + + allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS)) + + call phLR_A(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + call phLR_B(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + + call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) + + deallocate(Aph,Bph,XpY,XmY) ! Initialization @@ -87,48 +123,26 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH r_diis(:,:) = 0d0 rcond = 0d0 - amp(:,:,:) = 0d0 - res(:,:,:) = 0d0 + amp(:,:) = 0d0 + res(:,:) = 0d0 - ! Compute energy differences + ! Compute energy differences and coupling blocks - do i=nC+1,nO + do m=1,nS do j=nC+1,nO - do a=1,nV-nR - del(i,a,j) = eHF(i) + eHF(j) - eHF(nO+a) - eHF(p) + del(m,j) = Om(m) + eHF(j) - eHF(p) + vec(m,j) = sqrt(2d0)*rho(p,j,m) - end do end do end do - do i=nC+1,nO - do a=1,nV-nR - do b=1,nV-nR + do m=1,nS + do b=1,nV-nR - del(i,a,nO+b) = eHF(nO+a) + eHF(nO+b) - eHF(i) - eHF(p) + del(m,nO+b) = Om(m) + eHF(nO+b) - eHF(p) + vec(m,nO+b) = sqrt(2d0)*rho(p,nO+b,m) - end do - end do - end do - - do i=nC+1,nO - do a=1,nV-nR - do j=nC+1,nO - - vec(i,a,j) = sqrt(2d0)*ERI(p,nO+a,i,j) - - end do - end do - end do - - do i=nC+1,nO - do a=1,nV-nR - do b=1,nV-nR - - vec(i,a,nO+b) = sqrt(2d0)*ERI(p,i,nO+b,nO+a) - - end do end do end do @@ -152,41 +166,23 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH ! Compute residual for 2h1p sector - res(:,:,:) = vec(:,:,:) + (del(:,:,:) - Sig(p))*amp(:,:,:) + res(:,:) = vec(:,:) + (del(:,:) - Sig(p))*amp(:,:) - do i=nC+1,nO - do a=1,nV-nR - do j=nC+1,nO + do m=1,nS + do j=nC+1,nO - do k=nC+1,nO - do c=1,nV-nR + res(m,j) = res(m,j) - Om(m)*amp(m,j) - res(i,a,j) = res(i,a,j) - 2d0*ERI(j,nO+c,nO+a,k)*amp(i,c,k) & - - 2d0*ERI(i,nO+c,nO+a,k)*amp(k,c,j) - - end do - end do - - end do end do end do ! Compute residual for 2p1h sector - do i=nC+1,nO - do a=1,nV-nR - do b=1,nV-nR + do m=nC+1,nO + do b=1,nV-nR - do k=nC+1,nO - do c=1,nV-nR - - res(i,a,nO+b) = res(i,a,nO+b) + 2d0*ERI(nO+a,k,i,nO+c)*amp(k,c,nO+b) & - + 2d0*ERI(nO+b,k,i,nO+c)*amp(k,a,nO+c) - - end do - end do - - end do + res(m,nO+b) = res(m,nO+b) + Om(m)*amp(m,nO+b) + end do end do @@ -196,14 +192,14 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH ! Update amplitudes - amp(:,:,:) = amp(:,:,:) - res(:,:,:)/del(:,:,:) + amp(:,:) = amp(:,:) - res(:,:)/del(:,:) ! DIIS extrapolation if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nO*nV*nOrb,nO*nV*nOrb,n_diis,r_diis,t_diis,res,amp) + call DIIS_extrapolation(rcond,nS*nOrb,nS*nOrb,n_diis,r_diis,t_diis,res,amp) end if @@ -211,13 +207,11 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH Sig(p) = 0d0 - do i=nC+1,nO - do a=1,nV-nR - do q=nC+1,nOrb-nR + do m=1,nS + do q=nC+1,nOrb-nR - Sig(p) = Sig(p) + vec(i,a,q)*amp(i,a,q) + Sig(p) = Sig(p) + vec(m,q)*amp(m,q) - end do end do end do diff --git a/src/GW/ccRG0W0_TDA.f90 b/src/GW/ccRG0W0_TDA.f90 new file mode 100644 index 0000000..4c3ac4c --- /dev/null +++ b/src/GW/ccRG0W0_TDA.f90 @@ -0,0 +1,263 @@ +subroutine ccRG0W0_TDA(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + +! CC-based GW module (TDA screening) + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + double precision,intent(in) :: thresh + integer,intent(in) :: max_diis + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + 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) + +! Local variables + + integer :: p,q,r,s + integer :: i,j,k,l + integer :: a,b,c,d + + integer :: nSCF + double precision :: Conv + + double precision,allocatable :: Sig(:) + double precision,allocatable :: Z(:) + + double precision,allocatable :: del(:,:,:) + double precision,allocatable :: vec(:,:,:) + double precision,allocatable :: res(:,:,:) + double precision,allocatable :: amp(:,:,:) + + integer :: n_diis + double precision :: rcond + double precision,allocatable :: r_diis(:,:) + double precision,allocatable :: t_diis(:,:) + + +! Hello world + + write(*,*) + write(*,*)'*****************************' + write(*,*)'* CC-based G0W0 Calculation *' + write(*,*)'*****************************' + write(*,*) + +! Memory allocation + + allocate(del(nO,nV,nOrb)) + allocate(vec(nO,nV,nOrb)) + allocate(res(nO,nV,nOrb)) + allocate(amp(nO,nV,nOrb)) + + allocate(Sig(nOrb)) + allocate(Z(nOrb)) + + allocate(r_diis(nO*nV*nOrb,max_diis)) + allocate(t_diis(nO*nV*nOrb,max_diis)) + +! Initialization + + Sig(:) = 0d0 + Z(:) = 1d0 + +!-------------------------! +! Main loop over orbitals ! +!-------------------------! + + do p=nO,nO+1 + + ! Initialization + + Conv = 1d0 + nSCF = 0 + + n_diis = 0 + t_diis(:,:) = 0d0 + r_diis(:,:) = 0d0 + rcond = 0d0 + + amp(:,:,:) = 0d0 + res(:,:,:) = 0d0 + + ! Compute energy differences + + do i=nC+1,nO + do j=nC+1,nO + do a=1,nV-nR + + del(i,a,j) = eHF(i) + eHF(j) - eHF(nO+a) - eHF(p) + + end do + end do + end do + + do i=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + + del(i,a,nO+b) = eHF(nO+a) + eHF(nO+b) - eHF(i) - eHF(p) + + end do + end do + end do + + do i=nC+1,nO + do a=1,nV-nR + do j=nC+1,nO + + vec(i,a,j) = sqrt(2d0)*ERI(p,nO+a,i,j) + + end do + end do + end do + + do i=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + + vec(i,a,nO+b) = sqrt(2d0)*ERI(p,i,nO+b,nO+a) + + end do + end do + end do + + !----------------------! + ! Loop over amplitudes ! + !----------------------! + + write(*,*) + write(*,*)'-------------------------------------------------------------' + write(*,*)'| CC-based G0W0 calculation |' + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & + '|','#','|','Sig_c (eV)','|','e_GW (eV)','|','Conv','|' + write(*,*)'-------------------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Compute residual for 2h1p sector + + res(:,:,:) = vec(:,:,:) + (del(:,:,:) - Sig(p))*amp(:,:,:) + + do i=nC+1,nO + do a=1,nV-nR + do j=nC+1,nO + + do k=nC+1,nO + do c=1,nV-nR + + res(i,a,j) = res(i,a,j) - 2d0*ERI(j,nO+c,nO+a,k)*amp(i,c,k) +! - 2d0*ERI(i,nO+c,nO+a,k)*amp(k,c,j) + + end do + end do + + end do + end do + end do + + ! Compute residual for 2p1h sector + + do i=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + + do k=nC+1,nO + do c=1,nV-nR + + res(i,a,nO+b) = res(i,a,nO+b) + 2d0*ERI(nO+a,k,i,nO+c)*amp(k,c,nO+b) +! + 2d0*ERI(nO+b,k,i,nO+c)*amp(k,a,nO+c) + + end do + end do + + end do + end do + end do + + ! Check convergence + + Conv = maxval(abs(res)) + + ! Update amplitudes + + amp(:,:,:) = amp(:,:,:) - res(:,:,:)/del(:,:,:) + + ! DIIS extrapolation + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nV*nOrb,nO*nV*nOrb,n_diis,r_diis,t_diis,res,amp) + + end if + + ! Compute quasiparticle energy + + Sig(p) = 0d0 + + do i=nC+1,nO + do a=1,nV-nR + do q=nC+1,nOrb-nR + + Sig(p) = Sig(p) + vec(i,a,q)*amp(i,a,q) + + end do + end do + end do + + ! Dump results + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X)') & + '|',nSCF,'|',Sig(p)*HaToeV,'|',(eHF(p)+Sig(p))*HaToeV,'|',Conv,'|' + + end do + + write(*,*)'-------------------------------------------------------------' + write(*,*) + !------------------------------------------------------------------------ + ! End of SCF loop + !------------------------------------------------------------------------ + + ! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + end if + + write(*,*)'-------------------------------------------------------------------------------' + write(*,*)'| CC-based G0W0 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)') & + '|','Orb','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X)') & + '|',p,'|',eHF(p)*HaToeV,'|',Sig(p)*HaToeV,'|',Z(p),'|',(eHF(p)+Sig(p))*HaToeV,'|' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end do + +end subroutine diff --git a/src/GW/ccRG0W0_mat.f90 b/src/GW/ccRG0W0_mat.f90 deleted file mode 100644 index d75b68d..0000000 --- a/src/GW/ccRG0W0_mat.f90 +++ /dev/null @@ -1,278 +0,0 @@ -subroutine ccRG0W0_mat(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) - -! CC-based GW module - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: maxSCF - double precision,intent(in) :: thresh - - integer,intent(in) :: nBas - integer,intent(in) :: nOrb - integer,intent(in) :: nC - integer,intent(in) :: nO - integer,intent(in) :: nV - integer,intent(in) :: nR - 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) - -! Local variables - - integer :: p,q - integer :: i,j,k,l - integer :: a,b,c,d - - integer :: nSCF - double precision :: Conv - - double precision,allocatable :: delta_2h1p(:) - double precision,allocatable :: delta_2p1h(:) - - double precision,allocatable :: C_2h1p(:,:) - double precision,allocatable :: C_2p1h(:,:) - - double precision,allocatable :: V_2h1p(:) - double precision,allocatable :: V_2p1h(:) - - double precision,allocatable :: r_2h1p(:) - double precision,allocatable :: r_2p1h(:) - - double precision,allocatable :: t_2h1p(:) - double precision,allocatable :: t_2p1h(:) - - double precision,allocatable :: eGW(:) - double precision,allocatable :: Z(:) - - double precision,external :: Kronecker_delta - - integer :: n1h,n1p,n1h1p,n2h1p,n2p1h - integer :: ija,klc,iab,kcd - - -! Hello world - - write(*,*) - write(*,*)'*****************************' - write(*,*)'* CC-based G0W0 Calculation *' - write(*,*)'*****************************' - write(*,*) - -! Form energy denominator and guess amplitudes - - n1h = nO - n1p = nV - n1h1p = n1h*n1p - n2h1p = n1h*n1h1p - n2p1h = n1p*n1h1p - - allocate(delta_2h1p(n2h1p),delta_2p1h(n2p1h)) - allocate(C_2h1p(n2h1p,n2h1p),C_2p1h(n2p1h,n2p1h)) - allocate(V_2h1p(n2h1p),V_2p1h(n2p1h)) - allocate(t_2h1p(n2h1p),t_2p1h(n2p1h)) - allocate(r_2h1p(n2h1p),r_2p1h(n2p1h)) - allocate(eGW(nOrb),Z(nOrb)) - - ! Compute C2h1p and C2p1h - - ija = 0 - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nOrb-nR - ija = ija + 1 - - klc = 0 - do k=nC+1,nO - do l=nC+1,nO - do c=nO+1,nOrb-nR - klc = klc + 1 - - C_2h1p(ija,klc) = ((eHF(i) + eHF(j) - eHF(a))*Kronecker_delta(j,l)*Kronecker_delta(a,c) & - - 2d0*ERI(j,c,a,l))*Kronecker_delta(i,k) - - end do - end do - end do - - end do - end do - end do - - iab = 0 - do i=nC+1,nO - do a=nO+1,nOrb-nR - do b=nO+1,nOrb-nR - iab = iab + 1 - - kcd = 0 - do k=nC+1,nO - do c=nO+1,nOrb-nR - do d=nO+1,nOrb-nR - kcd = kcd + 1 - - C_2p1h(iab,kcd) = ((eHF(a) + eHF(b) - eHF(i))*Kronecker_delta(i,k)*Kronecker_delta(a,c) & - + 2d0*ERI(a,k,i,c))*Kronecker_delta(b,d) - - end do - end do - end do - - end do - end do - end do - -!-------------------------! -! Main loop over orbitals ! -!-------------------------! - - eGW(:) = eHF(:) - - do p=nO,nO - - ! Initialization - - Conv = 1d0 - nSCF = 0 - - t_2h1p(:) = 0d0 - t_2p1h(:) = 0d0 - - ! Compute energy differences - - ija = 0 - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nOrb-nR - ija = ija + 1 - - delta_2h1p(ija) = eHF(i) + eHF(j) - eHF(a) - eHF(p) - - end do - end do - end do - - iab = 0 - do i=nC+1,nO - do a=nO+1,nOrb-nR - do b=nO+1,nOrb-nR - iab = iab + 1 - - delta_2p1h(iab) = eHF(a) + eHF(b) - eHF(i) - eHF(p) - - end do - end do - end do - - klc = 0 - do k=nC+1,nO - do l=nC+1,nO - do c=nO+1,nOrb-nR - klc = klc + 1 - - V_2h1p(klc) = sqrt(2d0)*ERI(p,c,k,l) - - end do - end do - end do - - kcd = 0 - do k=nC+1,nO - do c=nO+1,nOrb-nR - do d=nO+1,nOrb-nR - kcd = kcd + 1 - - V_2p1h(kcd) = sqrt(2d0)*ERI(p,k,d,c) - - end do - end do - end do - - !----------------------! - ! Loop over amplitudes ! - !----------------------! - - write(*,*) - write(*,*)'----------------------------------------------' - write(*,*)'| CC-based G0W0 calculation |' - write(*,*)'----------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & - '|','#','|','HF','|','G0W0','|','Conv','|' - write(*,*)'----------------------------------------------' - - do while(Conv > thresh .and. nSCF < maxSCF) - - ! Increment - - nSCF = nSCF + 1 - - ! Compute residual for 2h1p sector - - r_2h1p = V_2h1p + matmul(C_2h1p,t_2h1p) - t_2h1p*eHF(p) & - - dot_product(t_2h1p,V_2h1p)*t_2h1p - t_2h1p*dot_product(V_2p1h,t_2p1h) - - ! Compute residual for 2p1h sector - - r_2p1h = V_2p1h + matmul(C_2p1h,t_2p1h) - t_2p1h*eHF(p) & - - t_2p1h*dot_product(V_2h1p,t_2h1p) - dot_product(t_2p1h,V_2p1h)*t_2p1h - - ! Check convergence - - Conv = max(maxval(abs(r_2h1p)),maxval(abs(r_2p1h))) - - ! Update amplitudes - - t_2h1p(:) = t_2h1p(:) - r_2h1p(:)/delta_2h1p(:) - t_2p1h(:) = t_2p1h(:) - r_2p1h(:)/delta_2p1h(:) - - ! Compute quasiparticle energies - - eGW(p) = eHF(p) + dot_product(V_2h1p,t_2h1p) + dot_product(V_2p1h,t_2p1h) - - ! Renormalization factor - - Z(:) = 1d0 - - ! Dump results - - write(*,'(1X,A1,1X,I3,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',eHF(p)*HaToeV,'|',eGW(p)*HaToeV,'|',Conv,'|' - - end do - - write(*,*)'----------------------------------------------' - !------------------------------------------------------------------------ - ! End of SCF loop - !------------------------------------------------------------------------ - - ! Did it actually converge? - - if(nSCF == maxSCF) then - - write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*) - - stop - - end if - - write(*,*)'-------------------------------------------------------------------------------' - write(*,*)' CC-G0W0 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)','|' - write(*,*)'-------------------------------------------------------------------------------' - - 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,'|',(eGW(p)-eHF(p))*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' - write(*,*)'-------------------------------------------------------------------------------' - - end do - -end subroutine