From 5d74b8073a3b3365181e584d48dd4a1d42f1b02e Mon Sep 17 00:00:00 2001 From: pfloos Date: Fri, 24 Nov 2023 15:31:29 +0100 Subject: [PATCH] cleaning up ufGW --- src/GW/GGW.f90 | 18 +- src/GW/RGW.f90 | 28 +-- src/GW/ufG0W0.f90 | 538 +++++++++++++++++++++++----------------------- src/GW/ufGW.f90 | 415 +++++++++++++++++++++++++---------- 4 files changed, 595 insertions(+), 404 deletions(-) diff --git a/src/GW/GGW.f90 b/src/GW/GGW.f90 index d6c6498..52e8ab7 100644 --- a/src/GW/GGW.f90 +++ b/src/GW/GGW.f90 @@ -1,6 +1,6 @@ -subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, & - TDA_W,TDA,dBSE,dTDA,linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & - ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF) +subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, & + TDA_W,TDA,dBSE,dTDA,linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc, & + ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,eHF) ! GW module @@ -45,8 +45,8 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan integer,intent(in) :: nR integer,intent(in) :: nS - double precision,intent(in) :: EHF - double precision,intent(in) :: epsHF(nBas2) + double precision,intent(in) :: EGHF + double precision,intent(in) :: eHF(nBas2) double precision,intent(in) :: cHF(nBas2,nBas2) double precision,intent(in) :: PHF(nBas2,nBas2) double precision,intent(in) :: S(nBas2,nBas2) @@ -71,7 +71,7 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan call wall_time(start_GW) call GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - linearize,eta,regularize,nBas2,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + linearize,eta,regularize,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -88,7 +88,7 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan call wall_time(start_GW) call evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - linearize,eta,regularize,nBas2,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + linearize,eta,regularize,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -105,8 +105,8 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan call wall_time(start_GW) call qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, & - dipole_int_AO,dipole_int,PHF,cHF,epsHF) + eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI, & + dipole_int_AO,dipole_int,PHF,cHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index cee59e7..f836401 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -1,7 +1,7 @@ -subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, & - exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & - linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & - ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF) +subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, & + exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & + linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & + ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,eHF) ! GW module @@ -50,8 +50,8 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre integer,intent(in) :: nR(nspin) integer,intent(in) :: nS(nspin) - double precision,intent(in) :: EHF - double precision,intent(in) :: epsHF(nBas) + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) @@ -76,7 +76,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre call wall_time(start_GW) call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -93,7 +93,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre call wall_time(start_GW) call evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -110,8 +110,8 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre call wall_time(start_GW) call qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, & - dipole_int_AO,dipole_int,PHF,cHF,epsHF) + singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI, & + dipole_int_AO,dipole_int,PHF,cHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -128,8 +128,8 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre call wall_time(start_GW) call SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA, & - singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, & - dipole_int_AO,dipole_int,PHF,cHF,epsHF) + singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI, & + dipole_int_AO,dipole_int,PHF,cHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -145,7 +145,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doufG0W0) then call wall_time(start_GW) - call ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF,TDA_W) + call ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -161,7 +161,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doufGW) then call wall_time(start_GW) - call ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF) + call ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/GW/ufG0W0.f90 b/src/GW/ufG0W0.f90 index 4bc1ea9..deaad60 100644 --- a/src/GW/ufG0W0.f90 +++ b/src/GW/ufG0W0.f90 @@ -1,4 +1,4 @@ -subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W) +subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Unfold G0W0 equations @@ -9,6 +9,7 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W) logical,intent(in) :: dotest + logical,intent(in) :: TDA_W integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -19,7 +20,6 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W) double precision,intent(in) :: ERHF double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: eHF(nBas) - logical,intent(in) :: TDA_W ! Local variables @@ -36,7 +36,6 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W) integer :: n2h1p,n2p1h,nH double precision,external :: Kronecker_delta double precision,allocatable :: H(:,:) - double precision,allocatable :: cGW(:,:) double precision,allocatable :: eGW(:) double precision,allocatable :: Z(:) double precision,allocatable :: Aph(:,:) @@ -55,9 +54,9 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W) ! Hello world write(*,*) - write(*,*)'**********************************************' - write(*,*)'| Unfolded G0W0 calculation |' - write(*,*)'**********************************************' + write(*,*)'****************************************' + write(*,*)'* Restricted Upfolded G0W0 Calculation *' + write(*,*)'****************************************' write(*,*) ! Dimension of the supermatrix @@ -68,7 +67,7 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W) ! Memory allocation - allocate(H(nH,nH),cGW(nH,nH),eGW(nH),Z(nH)) + allocate(H(nH,nH),eGW(nH),Z(nH)) ! Initialization @@ -77,230 +76,230 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W) H(:,:) = 0d0 - !!! Compute only the HOMO !!! - p=nO + p = nO if (TDA_W) then - ! TDA for W + ! TDA for W - write(*,*) 'Tamm-Dancoff approximation actived!' - write(*,*) + write(*,*) 'Tamm-Dancoff approximation actived!' + write(*,*) - !---------------------------! - ! Compute GW supermatrix ! - !---------------------------! - ! ! - ! | F V2h1p V2p1h | ! - ! | | ! - ! H = | V2h1p C2h1p 0 | ! - ! | | ! - ! | V2p1h 0 C2p1h | ! - ! ! - !---------------------------! + !---------------------------! + ! Compute GW supermatrix ! + !---------------------------! + ! ! + ! | F V2h1p V2p1h | ! + ! | | ! + ! H = | V2h1p C2h1p 0 | ! + ! | | ! + ! | V2p1h 0 C2p1h | ! + ! ! + !---------------------------! - !-------------! - ! Block C2h1p ! - !-------------! + !---------! + ! Block F ! + !---------! + + H(1,1) = eHF(p) - ija = 0 - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nBas-nR - ija = ija + 1 - - klc = 0 - do k=nC+1,nO - do l=nC+1,nO - do c=nO+1,nBas-nR - klc = klc + 1 - - H(1+ija,1+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 + !-------------! + ! Block V2h1p ! + !-------------! - end do - end do - end do - - !-------------! - ! Block C2p1h ! - !-------------! - - iab = 0 - do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - iab = iab + 1 - - kcd = 0 - do k=nC+1,nO - do c=nO+1,nBas-nR - do d=nO+1,nBas-nR - kcd = kcd + 1 - - H(1+n2h1p+iab,1+n2h1p+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 - - !---------! - ! Block F ! - !---------! - - H(1,1) = eHF(p) - - !-------------! - ! Block V2h1p ! - !-------------! - - klc = 0 - do k=nC+1,nO - do l=nC+1,nO - do c=nO+1,nBas-nR - klc = klc + 1 - - H(1 ,1+klc) = sqrt(2d0)*ERI(p,c,k,l) - H(1+klc,1 ) = sqrt(2d0)*ERI(p,c,k,l) - - end do - end do - end do - - !-------------! - ! Block V2p1h ! - !-------------! - - kcd = 0 - do k=nC+1,nO + klc = 0 + do k=nC+1,nO + do l=nC+1,nO do c=nO+1,nBas-nR - do d=nO+1,nBas-nR - - kcd = kcd + 1 - H(1 ,1+n2h1p+kcd) = sqrt(2d0)*ERI(p,k,d,c) - H(1+n2h1p+kcd,1 ) = sqrt(2d0)*ERI(p,k,d,c) - - end do + klc = klc + 1 + + H(1 ,1+klc) = sqrt(2d0)*ERI(p,c,k,l) + H(1+klc,1 ) = sqrt(2d0)*ERI(p,c,k,l) + end do - end do + end do + end do + + !-------------! + ! Block V2p1h ! + !-------------! + + kcd = 0 + do k=nC+1,nO + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + kcd = kcd + 1 + + H(1 ,1+n2h1p+kcd) = sqrt(2d0)*ERI(p,k,d,c) + H(1+n2h1p+kcd,1 ) = sqrt(2d0)*ERI(p,k,d,c) + + end do + end do + end do + + !-------------! + ! Block C2h1p ! + !-------------! + + ija = 0 + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + ija = ija + 1 + + klc = 0 + do k=nC+1,nO + do l=nC+1,nO + do c=nO+1,nBas-nR + klc = klc + 1 + + H(1+ija,1+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 + + !-------------! + ! Block C2p1h ! + !-------------! + + iab = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + iab = iab + 1 + + kcd = 0 + do k=nC+1,nO + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + kcd = kcd + 1 + + H(1+n2h1p+iab,1+n2h1p+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 else - ! RPA for W + ! RPA for W - write(*,*) 'Tamm-Dancoff approximation deactivated!' - write(*,*) + write(*,*) 'Tamm-Dancoff approximation deactivated!' + write(*,*) - !---------------------------! - ! Compute GW supermatrix ! - !---------------------------! - ! ! - ! | F W2h1p W2p1h | ! - ! | | ! - ! H = | W2h1p D2h1p 0 | ! - ! | | ! - ! | W2p1h 0 D2p1h | ! - ! ! - !---------------------------! + !---------------------------! + ! Compute GW supermatrix ! + !---------------------------! + ! ! + ! | F W2h1p W2p1h | ! + ! | | ! + ! H = | W2h1p D2h1p 0 | ! + ! | | ! + ! | W2p1h 0 D2p1h | ! + ! ! + !---------------------------! - ! Memory allocation ! - allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS)) + ! Memory allocation - ! Spin manifold + allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS)) - ispin = 1 + ! Spin manifold - !-------------------! - ! Compute screening ! - !-------------------! + ispin = 1 + + !-------------------! + ! Compute screening ! + !-------------------! call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) - call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) - !--------------------------! - ! Compute spectral weights ! - !--------------------------! + !--------------------------! + ! Compute spectral weights ! + !--------------------------! - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho) + call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho) - !---------! - ! Block F ! - !---------! + !---------! + ! Block F ! + !---------! - H(1,1) = eHF(p) + H(1,1) = eHF(p) - !-------------! - ! Block D2h1p ! - !-------------! + !-------------! + ! Block D2h1p ! + !-------------! - ija = 0 - do i=nC+1,nO - do ja=1,nS - ija = ija + 1 - H(1+ija,1+ija) = eHF(i) - Om(ja) - end do - end do + ija = 0 + do i=nC+1,nO + do ja=1,nS + ija = ija + 1 - !-------------! - ! Block W2h1p ! - !-------------! + H(1+ija,1+ija) = eHF(i) - Om(ja) - ija = 0 - do i=nC+1,nO - do ja=1,nS - ija = ija + 1 - H(1 ,1+ija) = sqrt(2d0)*rho(p,i,ja) - H(1+ija,1 ) = sqrt(2d0)*rho(p,i,ja) - end do - end do + end do + end do - !-------------! - ! Block D2h1p ! - !-------------! + !-------------! + ! Block W2h1p ! + !-------------! - iab = 0 - do b=nO+1,nBas-nR - ia = 0 - do i=nC+1,nO - do a=nO+1,nBas-nR - ia = ia + 1 - iab = iab + 1 - H(1+n2h1p+iab,1+n2h1p+iab) = eHF(b) + Om(ia) - end do - end do - end do + ija = 0 + do i=nC+1,nO + do ja=1,nS + ija = ija + 1 - !-------------! - ! Block W2p1h ! - !-------------! + H(1 ,1+ija) = sqrt(2d0)*rho(p,i,ja) + H(1+ija,1 ) = sqrt(2d0)*rho(p,i,ja) - iab = 0 - do b=nO+1,nBas-nR - ia = 0 - do i=nC+1,nO - do a=nO+1,nBas-nR - ia = ia + 1 - iab = iab + 1 - H(1 ,1+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia) - H(1+n2h1p+iab,1 ) = sqrt(2d0)*rho(p,b,ia) - end do - end do - end do + end do + end do + + !-------------! + ! Block D2p1h ! + !-------------! + + iab = 0 + do b=nO+1,nBas-nR + do ia=1,nS + iab = iab + 1 + + H(1+n2h1p+iab,1+n2h1p+iab) = eHF(b) + Om(ia) + + end do + end do + + !-------------! + ! Block W2p1h ! + !-------------! + + iab = 0 + do b=nO+1,nBas-nR + do ia=1,nS + iab = iab + 1 + + H(1 ,1+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia) + H(1+n2h1p+iab,1 ) = sqrt(2d0)*rho(p,b,ia) + + end do + end do end if @@ -308,96 +307,95 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W) ! Diagonalize supermatrix ! !-------------------------! - cGW(:,:) = H(:,:) - call diagonalize_matrix(nH,cGW,eGW) + H(:,:) = H(:,:) + call diagonalize_matrix(nH,H,eGW) !-----------------! ! Compute weights ! !-----------------! do s=1,nH - Z(s) = cGW(1,s)**2 + Z(s) = H(1,s)**2 end do !--------------! ! Dump results ! !--------------! - write(*,*)'-------------------------------------------' - write(*,'(A35,I3)')' G0W0 energies (eV) for orbital ',p - write(*,*)'-------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') & - '|','#','|','e_QP (eV)','|','Z','|' - write(*,*)'-------------------------------------------' + write(*,*)'-------------------------------------------' + write(*,'(1X,A32,I3,A8)')'| G0W0 energies (eV) for orbital',p,' |' + write(*,*)'-------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') & + '|','#','|','e_QP','|','Z','|' + write(*,*)'-------------------------------------------' - do s=1,nH - if(Z(s) > cutoff1) then - write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & - '|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|' - end if - end do - - write(*,*)'-------------------------------------------' - write(*,*) - - if(verbose) then - - do s=1,nH - - if(Z(s) > cutoff1) then - - write(*,*)'*************************************************************' - write(*,'(1X,A20,I3,A6,I3)')'Vector for orbital ',p,' and #',s - write(*,'(1X,A7,F10.6,A13,F10.6,1X)')' e_QP = ',eGW(s)*HaToeV,' eV and Z = ',Z(s) - write(*,*)'*************************************************************' - write(*,'(1X,A20,1X,A20,1X,A15,1X)') & - ' Configuration ',' Coefficient ',' Weight ' - write(*,*)'*************************************************************' - - if(p <= nO) & - write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') & - ' (',p,') ',cGW(1,s),cGW(1,s)**2 - if(p > nO) & - write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') & - ' (',p,') ',cGW(1,s),cGW(1,s)**2 - - klc = 0 - do k=nC+1,nO - do l=nC+1,nO - do c=nO+1,nBas-nR - - klc = klc + 1 - - if(abs(cGW(1+klc,s)) > cutoff2) & - write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') & - ' (',k,',',l,') -> (',c,') ',cGW(1+klc,s),cGW(1+klc,s)**2 - - end do - end do - end do - - kcd = 0 - do k=nC+1,nO - do c=nO+1,nBas-nR - do d=nO+1,nBas-nR - - kcd = kcd + 1 - if(abs(cGW(1+n2h1p+kcd,s)) > cutoff2) & - write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') & - ' (',k,') -> (',c,',',d,') ',cGW(1+n2h1p+kcd,s),cGW(1+n2h1p+kcd,s)**2 - - end do - end do - end do - - write(*,*)'*************************************************************' - write(*,*) - - end if - - end do - + do s=1,nH + if(Z(s) > cutoff1) then + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|' end if + end do + + write(*,*)'-------------------------------------------' + write(*,*) + if(verbose) then + + do s=1,nH + + if(Z(s) > cutoff1) then + + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A7,1X,I3,A6,I3,A1,1X,A7,F12.6,A13,F6.4,1X)') & + 'Orbital',p,' and #',s,':','e_QP = ',eGW(s)*HaToeV,' eV and Z = ',Z(s) + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A20,1X,A20,1X,A15,1X)') & + ' Configuration ',' Coefficient ',' Weight ' + write(*,*)'-------------------------------------------------------------' + + if(p <= nO) & + write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') & + ' (',p,') ',H(1,s),H(1,s)**2 + if(p > nO) & + write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') & + ' (',p,') ',H(1,s),H(1,s)**2 + + klc = 0 + do k=nC+1,nO + do l=nC+1,nO + do c=nO+1,nBas-nR + + klc = klc + 1 + + if(abs(H(1+klc,s)) > cutoff2) & + write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') & + ' (',k,',',l,') -> (',c,') ',H(1+klc,s),H(1+klc,s)**2 + + end do + end do + end do + + kcd = 0 + do k=nC+1,nO + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + kcd = kcd + 1 + if(abs(H(1+n2h1p+kcd,s)) > cutoff2) & + write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') & + ' (',k,') -> (',c,',',d,') ',H(1+n2h1p+kcd,s),H(1+n2h1p+kcd,s)**2 + + end do + end do + end do + + write(*,*)'-------------------------------------------------------------' + write(*,*) + + end if + + end do + + end if end subroutine diff --git a/src/GW/ufGW.f90 b/src/GW/ufGW.f90 index 66543ca..331b60b 100644 --- a/src/GW/ufGW.f90 +++ b/src/GW/ufGW.f90 @@ -1,4 +1,4 @@ -subroutine ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +subroutine ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Unfold GW equations @@ -9,6 +9,7 @@ subroutine ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) logical,intent(in) :: dotest + logical,intent(in) :: TDA_W integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -26,27 +27,36 @@ subroutine ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) integer :: s integer :: i,j,k,l integer :: a,b,c,d + integer :: ia,ja,kc,lc integer :: klc,kcd,ija,iab + logical :: dRPA + integer :: ispin + double precision :: EcRPA integer :: n2h1p,n2p1h,nH double precision,external :: Kronecker_delta double precision,allocatable :: H(:,:) double precision,allocatable :: eGW(:) double precision,allocatable :: Z(:) + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: rho(:,:,:) + + logical :: verbose = .true. + double precision,parameter :: cutoff1 = 0.01d0 + double precision,parameter :: cutoff2 = 0.01d0 ! Output variables ! Hello world write(*,*) - write(*,*)'**********************************************' - write(*,*)'| Unfolded GW calculation |' - write(*,*)'**********************************************' - write(*,*) - -! TDA for W - - write(*,*) 'Tamm-Dancoff approximation for dynamic screening by default!' + write(*,*)'**************************************' + write(*,*)'* Restricted Upfolded GW Calculation *' + write(*,*)'**************************************' write(*,*) ! Dimension of the supermatrix @@ -63,123 +73,248 @@ subroutine ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) H(:,:) = 0d0 -!---------------------------! -! Compute GW supermatrix ! -!---------------------------! -! ! -! | F V2h1p V2p1h | ! -! | | ! -! H = | V2h1p C2h1p 0 | ! -! | | ! -! | V2p1h 0 C2p1h | ! -! ! -!---------------------------! + if (TDA_W) then - !---------! - ! Block F ! - !---------! + ! TDA for W - do p=nC+1,nBas-nR - H(p,p) = eHF(p) - end do + write(*,*) 'Tamm-Dancoff approximation actived!' + write(*,*) - !-------------! - ! Block V2h1p ! - !-------------! + !---------------------------! + ! Compute GW supermatrix ! + !---------------------------! + ! ! + ! | F V2h1p V2p1h | ! + ! | | ! + ! H = | V2h1p C2h1p 0 | ! + ! | | ! + ! | V2p1h 0 C2p1h | ! + ! ! + !---------------------------! - do p=nC+1,nBas-nR - - klc = 0 - do k=nC+1,nO - do l=nC+1,nO + !---------! + ! Block F ! + !---------! + + do p=nC+1,nBas-nR + H(p,p) = eHF(p) + end do + + !-------------! + ! Block V2h1p ! + !-------------! + + do p=nC+1,nBas-nR + + klc = 0 + do k=nC+1,nO + do l=nC+1,nO + do c=nO+1,nBas-nR + klc = klc + 1 + + H(p ,nBas+klc) = sqrt(2d0)*ERI(p,c,k,l) + H(nBas+klc,p ) = sqrt(2d0)*ERI(p,c,k,l) + + end do + end do + end do + + end do + + !-------------! + ! Block V2p1h ! + !-------------! + + do p=nC+1,nBas-nR + + kcd = 0 + do k=nC+1,nO do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + kcd = kcd + 1 + + H(p ,nBas+n2h1p+kcd) = sqrt(2d0)*ERI(p,k,d,c) + H(nBas+n2h1p+kcd,p ) = sqrt(2d0)*ERI(p,k,d,c) + + end do + end do + end do + + end do + + !-------------! + ! Block C2h1p ! + !-------------! + + ija = 0 + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + ija = ija + 1 + + klc = 0 + do k=nC+1,nO + do l=nC+1,nO + do c=nO+1,nBas-nR + klc = klc + 1 + + H(nBas+ija,nBas+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 + + !-------------! + ! Block C2p1h ! + !-------------! + + iab = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + iab = iab + 1 + + kcd = 0 + do k=nC+1,nO + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + kcd = kcd + 1 + + H(nBas+n2h1p+iab,nBas+n2h1p+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 + + else + + ! RPA for W + + write(*,*) 'Tamm-Dancoff approximation deactivated!' + write(*,*) + + !---------------------------! + ! Compute GW supermatrix ! + !---------------------------! + ! ! + ! | F W2h1p W2p1h | ! + ! | | ! + ! H = | W2h1p D2h1p 0 | ! + ! | | ! + ! | W2p1h 0 D2p1h | ! + ! ! + !---------------------------! + + ! Memory allocation + + allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS)) + + ! Spin manifold + + ispin = 1 + + !-------------------! + ! Compute screening ! + !-------------------! + + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + + !--------------------------! + ! Compute spectral weights ! + !--------------------------! + + call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho) + + !---------! + ! Block F ! + !---------! + + do p=nC+1,nBas-nR + H(p,p) = eHF(p) + end do + + !-------------! + ! Block W2h1p ! + !-------------! + + do p=nC+1,nBas-nR + + klc = 0 + do k=nC+1,nO + do lc=1,nS klc = klc + 1 - H(p ,nBas+klc) = sqrt(2d0)*ERI(p,c,k,l) - H(nBas+klc,p ) = sqrt(2d0)*ERI(p,c,k,l) + H(p ,nBas+klc) = sqrt(2d0)*rho(p,k,lc) + H(nBas+klc,p ) = sqrt(2d0)*rho(p,k,lc) end do end do + end do - - end do - - !-------------! - ! Block V2p1h ! - !-------------! - - do p=nC+1,nBas-nR - - kcd = 0 - do k=nC+1,nO - do c=nO+1,nBas-nR + + !-------------! + ! Block W2p1h ! + !-------------! + + do p=nC+1,nBas-nR + + kcd = 0 + do kc=1,nS do d=nO+1,nBas-nR kcd = kcd + 1 - - H(p ,nBas+n2h1p+kcd) = sqrt(2d0)*ERI(p,k,d,c) - H(nBas+n2h1p+kcd,p ) = sqrt(2d0)*ERI(p,k,d,c) - + + H(p ,nBas+n2h1p+kcd) = sqrt(2d0)*rho(p,d,kc) + H(nBas+n2h1p+kcd,p ) = sqrt(2d0)*rho(p,d,kc) + end do end do + end do - - end do - - !-------------! - ! Block C2h1p ! - !-------------! - - ija = 0 - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nBas-nR + + !-------------! + ! Block D2h1p ! + !-------------! + + ija = 0 + do i=nC+1,nO + do ja=1,nS ija = ija + 1 - - klc = 0 - do k=nC+1,nO - do l=nC+1,nO - do c=nO+1,nBas-nR - klc = klc + 1 - - H(nBas+ija,nBas+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 - + + H(nBas+ija,nBas+ija) = eHF(i) - Om(ja) + end do end do - end do - - !-------------! - ! Block C2p1h ! - !-------------! - - iab = 0 - do i=nC+1,nO - do a=nO+1,nBas-nR + + !-------------! + ! Block D2p1h ! + !-------------! + + iab = 0 + do ia=1,nS do b=nO+1,nBas-nR iab = iab + 1 - - kcd = 0 - do k=nC+1,nO - do c=nO+1,nBas-nR - do d=nO+1,nBas-nR - kcd = kcd + 1 - - H(nBas+n2h1p+iab,nBas+n2h1p+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 - + + H(nBas+n2h1p+iab,nBas+n2h1p+iab) = eHF(b) + Om(ia) + end do end do - end do + + end if !-------------------------! ! Diagonalize supermatrix ! @@ -202,19 +337,77 @@ subroutine ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Dump results ! !--------------! - write(*,*)'-------------------------------------------' - write(*,*)' unfolded GW energies (eV) ' - write(*,*)'-------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') & - '|','#','|','e_QP (eV)','|','Z','|' - write(*,*)'-------------------------------------------' + write(*,*)'---------------------------------------------' + write(*,'(1X,A45)')'| GW energies (eV) for all orbitals |' + write(*,*)'---------------------------------------------' + write(*,'(1X,A1,1X,A5,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') & + '|','#','|','e_QP','|','Z','|' + write(*,*)'---------------------------------------------' do s=1,nH - write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + write(*,'(1X,A1,1X,I5,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|' enddo - write(*,*)'-------------------------------------------' + write(*,*)'---------------------------------------------' write(*,*) + if(verbose) then + + do s=1,nH + + if(Z(s) > cutoff1) then + + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A10,I5,A1,1X,A7,F12.6,A13,F6.4,1X)') & + 'Solution',s,':','e_QP = ',eGW(s)*HaToeV,' eV and Z = ',Z(s) + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A20,1X,A20,1X,A15,1X)') & + ' Configuration ',' Coefficient ',' Weight ' + write(*,*)'-------------------------------------------------------------' + + if(p <= nO) & + write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') & + ' (',p,') ',H(1,s),H(1,s)**2 + if(p > nO) & + write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') & + ' (',p,') ',H(1,s),H(1,s)**2 + + klc = 0 + do k=nC+1,nO + do l=nC+1,nO + do c=nO+1,nBas-nR + + klc = klc + 1 + + if(abs(H(1+klc,s)) > cutoff2) & + write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') & + ' (',k,',',l,') -> (',c,') ',H(1+klc,s),H(1+klc,s)**2 + + end do + end do + end do + + kcd = 0 + do k=nC+1,nO + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + kcd = kcd + 1 + if(abs(H(1+n2h1p+kcd,s)) > cutoff2) & + write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') & + ' (',k,') -> (',c,',',d,') ',H(1+n2h1p+kcd,s),H(1+n2h1p+kcd,s)**2 + + end do + end do + end do + + write(*,*)'-------------------------------------------------------------' + write(*,*) + + end if + + end do + + end if end subroutine