From ccd4ad570ba7c2915df85b5458e1a201e8cbe4de Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 26 Sep 2024 23:04:45 +0200 Subject: [PATCH] saving work in ppBSE@GW --- src/GW/GGW_ppBSE.f90 | 10 +- src/GW/GGW_ppBSE_dynamic_kernel_C.f90 | 34 +++-- src/GW/GGW_ppBSE_dynamic_perturbation.f90 | 2 - src/GW/GGW_ppBSE_static_kernel_C.f90 | 5 +- src/GW/GGW_ppBSE_upfolded.f90 | 162 ++++++++++++++++------ 5 files changed, 143 insertions(+), 70 deletions(-) diff --git a/src/GW/GGW_ppBSE.f90 b/src/GW/GGW_ppBSE.f90 index e6f0d1a..5c48782 100644 --- a/src/GW/GGW_ppBSE.f90 +++ b/src/GW/GGW_ppBSE.f90 @@ -65,22 +65,18 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int, ! Compute RPA screening ! !-----------------------! - EcRPA = 0d0 - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), & Aph(nS,nS),Bph(nS,nS)) call phGLR_A(dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph) if(.not.TDA_W) call phGLR_B(dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) - call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phGLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call GGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) deallocate(XpY_RPA,XmY_RPA,Aph,Bph) - EcBSE = 0d0 - nOO = nO*(nO-1)/2 nVV = nV*(nV-1)/2 @@ -103,7 +99,7 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int, Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) - call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE) + call ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE) call ppLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) diff --git a/src/GW/GGW_ppBSE_dynamic_kernel_C.f90 b/src/GW/GGW_ppBSE_dynamic_kernel_C.f90 index 76f02ce..75c0050 100644 --- a/src/GW/GGW_ppBSE_dynamic_kernel_C.f90 +++ b/src/GW/GGW_ppBSE_dynamic_kernel_C.f90 @@ -53,30 +53,34 @@ subroutine GGW_ppBSE_dynamic_kernel_C(eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eGW,Om, do m=1,nS - dem = OmBSE - eGW(c) - Om(m) - eGW(b) + dem = OmBSE - (eGW(a) + eGW(c) + Om(m)) +! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m)) + num = - rho(b,c,m)*rho(a,d,m) + + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + dem = OmBSE - (eGW(b) + eGW(d) + Om(m)) +! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m)) + num = - rho(b,c,m)*rho(a,d,m) + + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + dem = OmBSE - (eGW(b) + eGW(c) + Om(m)) num = rho(a,c,m)*rho(b,d,m) +! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m)) KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - dem = OmBSE - eGW(c) - Om(m) - eGW(a) - num = rho(b,c,m)*rho(a,d,m) - - KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - - dem = OmBSE - eGW(d) - Om(m) - eGW(a) - num = rho(a,c,m)*rho(b,d,m) + dem = OmBSE - (eGW(a) + eGW(d) + Om(m)) + num = rho(a,c,m)*rho(b,d,m) +! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m)) KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - dem = OmBSE - eGW(d) - Om(m) - eGW(b) - num = rho(b,c,m)*rho(a,d,m) - - KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - end do end do diff --git a/src/GW/GGW_ppBSE_dynamic_perturbation.f90 b/src/GW/GGW_ppBSE_dynamic_perturbation.f90 index a4113e6..344c207 100644 --- a/src/GW/GGW_ppBSE_dynamic_perturbation.f90 +++ b/src/GW/GGW_ppBSE_dynamic_perturbation.f90 @@ -71,8 +71,6 @@ subroutine GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,e write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------' - print*,nVV,maxVV - do ab=1,min(nVV,maxVV) if(dTDA) then diff --git a/src/GW/GGW_ppBSE_static_kernel_C.f90 b/src/GW/GGW_ppBSE_static_kernel_C.f90 index 0e8cd05..259951c 100644 --- a/src/GW/GGW_ppBSE_static_kernel_C.f90 +++ b/src/GW/GGW_ppBSE_static_kernel_C.f90 @@ -37,14 +37,11 @@ subroutine GGW_ppBSE_static_kernel_C(eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI,Om,r double precision,intent(out) :: KC(nVV,nVV) -!---------------! -! SpinOrb block ! -!---------------! - ab = 0 do a=nO+1,nBas-nR do b=a+1,nBas-nR ab = ab + 1 + cd = 0 do c=nO+1,nBas-nR do d=c+1,nBas-nR diff --git a/src/GW/GGW_ppBSE_upfolded.f90 b/src/GW/GGW_ppBSE_upfolded.f90 index 928d175..85d757c 100644 --- a/src/GW/GGW_ppBSE_upfolded.f90 +++ b/src/GW/GGW_ppBSE_upfolded.f90 @@ -24,7 +24,8 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) integer :: i,j,k,l integer :: a,b,c,d integer :: m,ij,kl,ijm - integer,parameter :: maxH = 1000 + integer :: ab,cd,abm + integer,parameter :: maxH = 100 double precision :: tmp,tmp1,tmp2,tmp3,tmp4 integer :: n1h,n1p,n2h,n2p,n1h1p,n3h1p,n3p1h,n2h2p,nH @@ -64,7 +65,8 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) n3h1p = n2h*n1h1p n3p1h = n2p*n1h1p - nH = n2h + n3h1p + nH = n2p + n3p1h +! nH = n2h + n3h1p ! Memory allocation @@ -91,21 +93,21 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) !----------------------------------------! !---------! - ! Block D ! + ! Block C ! !---------! - ij = 0 - do i=nC+1,nO - do j=i+1,nO - ij = ij + 1 + ab = 0 + do a=nO+1,nOrb-nR + do b=a+1,nOrb-nR + ab = ab + 1 - kl = 0 - do k=nC+1,nO - do l=k+1,nO - kl = kl + 1 + cd = 0 + do c=nO+1,nOrb-nR + do d=c+1,nOrb-nR + cd = cd + 1 - H(ij,kl) = - (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - + (ERI(i,j,k,l) - ERI(i,j,l,k)) + H(ab,cd) = (eGW(a) + eGW(b))*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + + (ERI(a,b,c,d) - ERI(a,b,d,c)) end do end do @@ -113,29 +115,52 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) end do end do - !----------------! - ! Blocks M1 & M2 ! - !----------------! + !---------! + ! Block D ! + !---------! - ijm = 0 - do i=nC+1,nO - do j=i+1,nO +! ij = 0 +! do i=nC+1,nO +! do j=i+1,nO +! ij = ij + 1 +! +! kl = 0 +! do k=nC+1,nO +! do l=k+1,nO +! kl = kl + 1 +! +! H(ij,kl) = - (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) & +! + (ERI(i,j,k,l) - ERI(i,j,l,k)) +! +! end do +! end do +! +! end do +! end do + + !-----------------! + ! Coupling Blocks ! + !-----------------! + + abm = 0 + do a=nO+1,nOrb-nR + do b=a+1,nOrb-nR do m=1,nS - ijm = ijm + 1 + abm = abm + 1 - kl = 0 - do k=nC+1,nO - do l=k+1,nO - kl = kl + 1 + cd = 0 + do c=nO+1,nOrb-nR + do d=c+1,nOrb-nR + cd = cd + 1 - tmp1 = Kronecker_delta(j,l)*rho(i,k,m) - tmp2 = Kronecker_delta(j,k)*rho(i,l,m) - tmp3 = Kronecker_delta(i,l)*rho(j,k,m) - tmp4 = Kronecker_delta(i,k)*rho(j,l,m) + tmp1 = Kronecker_delta(b,d)*rho(a,c,m) + tmp2 = Kronecker_delta(b,c)*rho(a,d,m) + tmp3 = Kronecker_delta(a,d)*rho(b,c,m) + tmp4 = Kronecker_delta(a,c)*rho(b,d,m) + + H(n2p+0*n3p1h+abm,cd ) = tmp1 + tmp2 + H(cd ,n2p+0*n3p1h+abm) = tmp3 + tmp4 - H(n2h+0*n3h1p+ijm,kl ) = tmp1 - tmp2 - H(kl ,n2h+0*n3h1p+ijm) = tmp3 - tmp4 - ! H(n2h+1*n3h1p+ijm,kl ) = +tmp4 ! H(kl ,n2h+1*n3h1p+ijm) = +tmp2 ! @@ -152,27 +177,80 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) end do end do +! ijm = 0 +! do i=nC+1,nO +! do j=i+1,nO +! do m=1,nS +! ijm = ijm + 1 + +! kl = 0 +! do k=nC+1,nO +! do l=k+1,nO +! kl = kl + 1 + +! tmp1 = Kronecker_delta(j,l)*rho(i,k,m) +! tmp2 = Kronecker_delta(j,k)*rho(i,l,m) +! tmp3 = Kronecker_delta(i,l)*rho(j,k,m) +! tmp4 = Kronecker_delta(i,k)*rho(j,l,m) + +! H(n2h+0*n3h1p+ijm,kl ) = tmp1 - tmp2 +! H(kl ,n2h+0*n3h1p+ijm) = tmp3 - tmp4 +! +! H(n2h+1*n3h1p+ijm,kl ) = +tmp4 +! H(kl ,n2h+1*n3h1p+ijm) = +tmp2 +! +! H(n2h+2*n3h1p+ijm,kl ) = +tmp1 +! H(kl ,n2h+2*n3h1p+ijm) = +tmp4 + +! H(n2h+3*n3h1p+ijm,kl ) = +tmp3 +! H(kl ,n2h+3*n3h1p+ijm) = +tmp1 + +! end do +! end do + +! end do +! end do +! end do + !------------! - ! Block 3h1p ! + ! Block 3p1h ! !------------! - ijm = 0 - do i=nC+1,nO - do j=i+1,nO + abm = 0 + do a=nO+1,nOrb-nR + do b=a+1,nOrb-nR do m=1,nS - ijm = ijm + 1 + abm = abm + 1 - tmp = - eGW(i) - eGW(j) + Om(m) + tmp = eGW(a) + eGW(b) + Om(m) - H(n2h+0*n3h1p+ijm,n2h+0*n3h1p+ijm) = tmp -! H(n2h+1*n3h1p+ijm,n2h+1*n3h1p+ijm) = tmp -! H(n2h+2*n3h1p+ijm,n2h+2*n3h1p+ijm) = tmp -! H(n2h+3*n3h1p+ijm,n2h+3*n3h1p+ijm) = tmp + H(n2p+0*n3p1h+abm,n2p+0*n3p1h+abm) = tmp end do end do end do + !------------! + ! Block 3h1p ! + !------------! + +! ijm = 0 +! do i=nC+1,nO +! do j=i+1,nO +! do m=1,nS +! ijm = ijm + 1 + +! tmp = - eGW(i) - eGW(j) + Om(m) + +! H(n2h+0*n3h1p+ijm,n2h+0*n3h1p+ijm) = tmp +! H(n2h+1*n3h1p+ijm,n2h+1*n3h1p+ijm) = tmp +! H(n2h+2*n3h1p+ijm,n2h+2*n3h1p+ijm) = tmp +! H(n2h+3*n3h1p+ijm,n2h+3*n3h1p+ijm) = tmp + +! end do +! end do +! end do + !-------------------------! ! Diagonalize supermatrix ! !-------------------------! @@ -211,7 +289,7 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) 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,'|',-OmBSE(s)*HaToeV,'|',Z(s),'|' + '|',s,'|',OmBSE(s)*HaToeV,'|',Z(s),'|' end do write(*,*)'-------------------------------------------'