diff --git a/src/GT/GGTpp_self_energy_diag.f90 b/src/GT/GGTpp_self_energy_diag.f90 index f319548..961c915 100644 --- a/src/GT/GGTpp_self_energy_diag.f90 +++ b/src/GT/GGTpp_self_energy_diag.f90 @@ -43,129 +43,129 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh ! Occupied part of the Tpp self-energy ! !--------------------------------------! - do p=nC+1,nBas-nR - do i=nC+1,nO +! do p=nC+1,nBas-nR +! do i=nC+1,nO - do cd=1,nVV - eps = e(p) + e(i) - Om1(cd) - num = rho1(p,i,cd)**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 +! do cd=1,nVV +! eps = e(p) + e(i) - Om1(cd) +! num = rho1(p,i,cd)**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 - end do +! end do +! end do -!------------------------------------------! -! Virtual part of the T-matrix self-energy ! -!------------------------------------------! +! !------------------------------------------! +! ! Virtual part of the T-matrix self-energy ! +! !------------------------------------------! - do p=nC+1,nBas-nR - do a=nO+1,nBas-nR +! do p=nC+1,nBas-nR +! do a=nO+1,nBas-nR - do kl=1,nOO - eps = e(p) + e(a) - Om2(kl) - num = rho2(p,a,kl)**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 +! do kl=1,nOO +! eps = e(p) + e(a) - Om2(kl) +! num = rho2(p,a,kl)**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 - end do +! end do +! end do !-----------------------------------------------! ! Testing another way to compute GT self-energy ! !-----------------------------------------------! - ! do p=nC+1,nBas-nR - ! do i=nC+1,nO - ! do j=nC+1,nO - ! do a=nO+1,nBas-nR + do p=nC+1,nBas-nR + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR - ! eps = e(p) + e(a) - e(i) - e(j) - ! num = 0.5d0*(ERI(p,a,i,j) - ERI(p,a,j,i))**2 + eps = e(p) + e(a) - e(i) - e(j) + num = 0.5d0*(ERI(p,a,i,j) - ERI(p,a,j,i))**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 + 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 - ! do a=nO+1,nBas-nR + end do + do a=nO+1,nBas-nR - ! do m=1,nVV - ! num = - ERI(p,a,i,j) * rho1(p,a,m) * rho1(i,j,m) - ! dem1 = e(p) + e(a) - e(i) - e(j) - ! dem2 = Om1(m) - e(i) - e(j) - ! Sig(p) = Sig(p) + num/dem1/dem2 - ! Z(p) = Z(p) - num/dem1/dem1/dem2 - ! end do + do m=1,nVV + num = - ERI(p,a,i,j) * rho1(p,a,m) * rho1(i,j,m) + dem1 = e(p) + e(a) - e(i) - e(j) + dem2 = Om1(m) - e(i) - e(j) + Sig(p) = Sig(p) + num/dem1/dem2 + Z(p) = Z(p) - num/dem1/dem1/dem2 + end do - ! do m=1,nOO - ! num = - ERI(p,a,i,j) * rho2(p,a,m) * rho2(i,j,m) - ! dem1 = e(p) + e(a) - e(i) - e(j) - ! dem2 = e(p) + e(a) - Om2(m) - ! Sig(p) = Sig(p) + num/dem1/dem2 - ! Z(p) = Z(p) - num/dem1/dem1/dem2 - num/dem1/dem2/dem2 - ! end do + do m=1,nOO + num = - ERI(p,a,i,j) * rho2(p,a,m) * rho2(i,j,m) + dem1 = e(p) + e(a) - e(i) - e(j) + dem2 = e(p) + e(a) - Om2(m) + Sig(p) = Sig(p) + num/dem1/dem2 + Z(p) = Z(p) - num/dem1/dem1/dem2 - num/dem1/dem2/dem2 + end do - ! end do - ! do k=nC+1,nO + end do + ! do k=nC+1,nO - ! do m=1,nVV - ! num = - ERI(p,i,j,k) * rho1(p,i,m) * rho1(j,k,m) - ! dem1 = e(p) + e(i) - Om1(m) - ! dem2 = Om1(m) - e(j) - e(k) - ! Sig(p) = Sig(p) + num/dem1/dem2 - ! Z(p) = Z(p) - num/dem1/dem1/dem2 - ! end do + ! do m=1,nVV + ! num = - ERI(p,i,j,k) * rho1(p,i,m) * rho1(j,k,m) + ! dem1 = e(p) + e(i) - Om1(m) + ! dem2 = Om1(m) - e(j) - e(k) + ! Sig(p) = Sig(p) + num/dem1/dem2 + ! Z(p) = Z(p) - num/dem1/dem1/dem2 + ! end do - ! end do - ! end do - ! end do - ! end do - ! do p=nC+1,nBas-nR - ! do a=nO+1,nBas-nR - ! do b=nO+1,nBas-nR - ! do i=nC+1,nO + ! end do + end do + end do + end do + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + do i=nC+1,nO - ! eps = e(p) + e(i) - e(a) - e(b) - ! num = 0.5d0*(ERI(p,i,a,b) - ERI(p,i,b,a))**2 + eps = e(p) + e(i) - e(a) - e(b) + num = 0.5d0*(ERI(p,i,a,b) - ERI(p,i,b,a))**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 + 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 - ! do i=nC+1,nO + end do + do i=nC+1,nO - ! do m=1,nVV - ! num = ERI(p,i,a,b) * rho1(p,i,m) * rho1(a,b,m) - ! dem1 = e(p) + e(i) - e(a) - e(b) - ! dem2 = e(p) + e(i) - Om1(m) - ! Sig(p) = Sig(p) + num/dem1/dem2 - ! Z(p) = Z(p) - num/dem1/dem1/dem2 - num/dem1/dem2/dem2 - ! end do + do m=1,nVV + num = ERI(p,i,a,b) * rho1(p,i,m) * rho1(a,b,m) + dem1 = e(p) + e(i) - e(a) - e(b) + dem2 = e(p) + e(i) - Om1(m) + Sig(p) = Sig(p) + num/dem1/dem2 + Z(p) = Z(p) - num/dem1/dem1/dem2 - num/dem1/dem2/dem2 + end do - ! do m=1,nOO - ! num = ERI(p,i,a,b) * rho2(p,i,m) * rho2(a,b,m) - ! dem1 = e(p) + e(i) - e(a) - e(b) - ! dem2 = Om2(m) - e(a) - e(b) - ! Sig(p) = Sig(p) + num/dem1/dem2 - ! Z(p) = Z(p) - num/dem1/dem1/dem2 - ! end do + do m=1,nOO + num = ERI(p,i,a,b) * rho2(p,i,m) * rho2(a,b,m) + dem1 = e(p) + e(i) - e(a) - e(b) + dem2 = Om2(m) - e(a) - e(b) + Sig(p) = Sig(p) + num/dem1/dem2 + Z(p) = Z(p) - num/dem1/dem1/dem2 + end do - ! end do - ! do c=nO+1,nBas-nR - ! do m=1,nOO - ! num = ERI(p,a,b,c) * rho2(p,a,m) * rho2(b,c,m) - ! dem1 = e(p) + e(a) - Om2(m) - ! dem2 = Om2(m) - e(b) - e(c) - ! Sig(p) = Sig(p) + num/dem1/dem2 - ! Z(p) = Z(p) - num/dem1/dem1/dem2 - ! end do + end do + ! do c=nO+1,nBas-nR + ! do m=1,nOO + ! num = ERI(p,a,b,c) * rho2(p,a,m) * rho2(b,c,m) + ! dem1 = e(p) + e(a) - Om2(m) + ! dem2 = Om2(m) - e(b) - e(c) + ! Sig(p) = Sig(p) + num/dem1/dem2 + ! Z(p) = Z(p) - num/dem1/dem1/dem2 + ! end do - ! end do - ! end do - ! end do - ! end do + ! end do + end do + end do + end do !-------------------------------------! ! Galitskii-Migdal correlation energy ! diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90 index a58989c..72163d6 100644 --- a/src/GW/GG0W0.f90 +++ b/src/GW/GG0W0.f90 @@ -119,7 +119,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA if(doSRG) then call GGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) else - call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) + call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z,ERI) end if !-----------------------------------! diff --git a/src/GW/GGW_self_energy_diag.f90 b/src/GW/GGW_self_energy_diag.f90 index 2a02837..0f1bd30 100644 --- a/src/GW/GGW_self_energy_diag.f90 +++ b/src/GW/GGW_self_energy_diag.f90 @@ -1,4 +1,4 @@ -subroutine GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) +subroutine GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z,ERI) ! Compute diagonal of the correlation part of the self-energy and the renormalization factor @@ -17,11 +17,12 @@ subroutine GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) double precision,intent(in) :: e(nBas) double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables - integer :: i,a,p,m - double precision :: num,eps + integer :: i,j,a,b,p,m + double precision :: num,eps,dem1,dem2 ! Output variables @@ -38,36 +39,118 @@ subroutine GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) ! GW self-energy ! !----------------! -! Occupied part of the correlation self-energy +! ! Occupied part of the correlation self-energy + +! do p=nC+1,nBas-nR +! do i=nC+1,nO +! do m=1,nS + +! eps = e(p) - e(i) + Om(m) +! num = rho(p,i,m)**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 +! end do + +! ! Virtual part of the correlation self-energy + +! do p=nC+1,nBas-nR +! do a=nO+1,nBas-nR +! do m=1,nS + +! eps = e(p) - e(a) - Om(m) +! num = rho(p,a,m)**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 +! end do + +!-----------------------------------------------! +! Testing another way to compute GT self-energy ! +!-----------------------------------------------! do p=nC+1,nBas-nR - do i=nC+1,nO - do m=1,nS + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR - eps = e(p) - e(i) + Om(m) - num = rho(p,i,m)**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 + eps = e(p) + e(a) - e(i) - e(j) + num = ERI(p,a,i,j)**2 - end do - end do + 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 + do a=nO+1,nBas-nR + + do m=1,nS + num = - ERI(p,i,j,a) * rho(i,a,m) * rho(j,p,m) + dem1 = e(p) + e(a) - e(i) - e(j) + dem2 = e(p) - e(j) + Om(m) + Sig(p) = Sig(p) + num/dem1/dem2 + Z(p) = Z(p) - num/dem1/dem1/dem2 - num/dem1/dem2/dem2 + + num = - ERI(p,i,j,a) * rho(i,a,m) * rho(j,p,m) + dem1 = e(p) + e(a) - e(i) - e(j) + dem2 = e(a) - e(i) + Om(m) + Sig(p) = Sig(p) + num/dem1/dem2 + Z(p) = Z(p) - num/dem1/dem1/dem2 + + num = - ERI(p,a,j,i) * rho(a,i,m) * rho(j,p,m) + dem1 = e(p) - e(j) + Om(m) + dem2 = e(a) - e(i) + Om(m) + Sig(p) = Sig(p) + num/dem1/dem2 + Z(p) = Z(p) - num/dem1/dem1/dem2 + end do + + end do + end do + end do end do - -! Virtual part of the correlation self-energy - do p=nC+1,nBas-nR - do a=nO+1,nBas-nR - do m=1,nS + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + do i=nC+1,nO + + eps = e(p) + e(i) - e(a) - e(b) + num = ERI(p,i,a,b)**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 + do i=nC+1,nO + + do m=1,nS + num = ERI(p,a,b,i) * rho(a,i,m) * rho(b,p,m) + dem1 = e(p) + e(i) - e(a) - e(b) + dem2 = e(p) - e(b) - Om(m) + Sig(p) = Sig(p) + num/dem1/dem2 + Z(p) = Z(p) - num/dem1/dem1/dem2 - num/dem1/dem2/dem2 - eps = e(p) - e(a) - Om(m) - num = rho(p,a,m)**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 + num = - ERI(p,a,b,i) * rho(a,i,m) * rho(b,p,m) + dem1 = e(p) + e(i) - e(a) - e(b) + dem2 = e(a) - e(i) + Om(m) + Sig(p) = Sig(p) + num/dem1/dem2 + Z(p) = Z(p) - num/dem1/dem1/dem2 + + num = - ERI(p,i,b,a) * rho(i,a,m) * rho(b,p,m) + dem1 = e(p) - e(b) - Om(m) + dem2 = e(a) - e(i) + Om(m) + Sig(p) = Sig(p) + num/dem1/dem2 + Z(p) = Z(p) - num/dem1/dem1/dem2 + end do - end do - end do + end do + end do + end do end do + ! Galitskii-Migdal correlation energy EcGM = 0d0 diff --git a/src/Parquet/G_Parquet_self_energy.f90 b/src/Parquet/G_Parquet_self_energy.f90 index 7ab19b5..d163a42 100644 --- a/src/Parquet/G_Parquet_self_energy.f90 +++ b/src/Parquet/G_Parquet_self_energy.f90 @@ -55,6 +55,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& eps = eQP(p) + eQP(a) - eQP(i) - eQP(j) reg = (1d0 - exp(- 2d0 * eta * eps * eps)) num = 0.5d0*(ERI(p,a,j,i) - ERI(p,a,i,j))**2 + ! num = ERI(p,a,j,i)**2 SigC(p) = SigC(p) + num*reg/eps Z(p) = Z(p) - num*reg/eps**2 @@ -70,6 +71,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& eps = eQP(p) + eQP(i) - eQP(a) - eQP(b) reg = (1d0 - exp(- 2d0 * eta * eps * eps)) num = 0.5d0*(ERI(p,i,b,a) - ERI(p,i,a,b))**2 + ! num = ERI(p,i,b,a)**2 SigC(p) = SigC(p) + num*reg/eps Z(p) = Z(p) - num*reg/eps**2 @@ -89,151 +91,108 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& ! eh part of the self-energy ! !-----------------------------! - ! call wall_time(start_t) + call wall_time(start_t) - ! !$OMP PARALLEL DEFAULT(NONE) & - ! !$OMP PRIVATE(p,i,a,j,b,n,num,dem1,dem2,reg1,reg2) & - ! !$OMP SHARED(nC,nO,nOrb,nR,nS,eta,ERI,eQP,eh_rho,eh_Om,SigC,Z) - ! !$OMP DO COLLAPSE(2) - ! do p=nC+1,nOrb-nR + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(p,i,a,j,b,n,num,dem1,dem2,reg1,reg2) & + !$OMP SHARED(nC,nO,nOrb,nR,nS,eta,ERI,eQP,eh_rho,eh_Om,SigC,Z) + !$OMP DO COLLAPSE(2) + do p=nC+1,nOrb-nR - ! do i=nC+1,nO - ! do a=nO+1,nOrb-nR + do i=nC+1,nO + do a=nO+1,nOrb-nR - ! do n=1,nS - ! !3h2p - ! do j=nC+1,nO + do n=1,nS + !3h2p + do j=nC+1,nO - ! num = - ERI(p,a,j,i) * eh_rho(p,j,nS+n) * eh_rho(i,a,nS+n) & - ! + ERI(p,a,i,j) * eh_rho(a,i,n) * eh_rho(j,p,n) - ! dem1 = eQP(p) - eQP(j) + eh_Om(n) - ! dem2 = eQP(p) - eQP(j) + eQP(a) - eQP(i) - ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + num = - (ERI(p,a,j,i) - ERI(p,a,i,j)) * eh_rho(p,j,nS+n) * eh_rho(i,a,nS+n) !& + !+ ERI(p,a,i,j) * eh_rho(a,i,n) * eh_rho(j,p,n) + dem1 = eQP(p) - eQP(j) + eh_Om(n) + dem2 = eQP(p) - eQP(j) + eQP(a) - eQP(i) + reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) & - ! - num * (reg1/dem1/dem1) * (reg2/dem2) + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) & + - num * (reg1/dem1/dem1) * (reg2/dem2) - ! ! num = - 0d0*ERI(p,a,j,i) * eh_rho(p,j,nS+n) * eh_rho(i,a,nS+n) & - ! ! + ERI(p,a,i,j) * eh_rho(a,i,n) * eh_rho(j,p,n) - ! ! dem1 = eQP(a) - eQP(i) - eh_Om(n) - ! ! dem2 = eQP(p) - eQP(j) + eh_Om(n) - ! ! reg1 = (1d0 - 0d0*exp(- 2d0 * eta * dem1 * dem1)) - ! ! reg2 = (1d0 - 0d0*exp(- 2d0 * eta * dem2 * dem2)) + num = - (ERI(p,i,j,a) - ERI(p,i,a,j)) * eh_rho(p,j,nS+n) * eh_rho(a,i,nS+n) !& + !+ ERI(p,i,a,j) * eh_rho(i,a,n) * eh_rho(j,p,n) - ! ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - - ! ! num = + ERI(p,a,j,i) * eh_rho(p,j,nS+n) * eh_rho(i,a,nS+n) & - ! ! - ERI(p,a,i,j) * eh_rho(a,i,n) * eh_rho(j,p,n) - - ! ! dem1 = eQP(a) - eQP(i) - eh_Om(n) - ! ! dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j) - ! ! reg1 = (1d0 - 0d0*exp(- 2d0 * eta * dem1 * dem1)) - ! ! reg2 = (1d0 - 0d0*exp(- 2d0 * eta * dem2 * dem2)) - - ! ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - - ! num = - ERI(p,i,j,a) * eh_rho(p,j,nS+n) * eh_rho(a,i,nS+n) & - ! + ERI(p,i,a,j) * eh_rho(i,a,n) * eh_rho(j,p,n) - - ! dem1 = eQP(a) - eQP(i) + eh_Om(n) - ! dem2 = eQP(p) - eQP(j) + eh_Om(n) - ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - - ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) + dem1 = eQP(a) - eQP(i) + eh_Om(n) + dem2 = eQP(p) - eQP(j) + eh_Om(n) + reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - ! num = - ERI(p,a,j,i) * eh_rho(p,j,n) * eh_rho(i,a,n) & - ! + ERI(p,a,i,j) * eh_rho(a,i,nS+n) * eh_rho(j,p,nS+n) + num = - (ERI(p,a,j,i) - ERI(p,a,i,j)) * eh_rho(p,j,n) * eh_rho(i,a,n) !& + !+ ERI(p,a,i,j) * eh_rho(a,i,nS+n) * eh_rho(j,p,nS+n) - ! dem1 = eQP(a) - eQP(i) + eh_Om(n) - ! dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j) - ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + dem1 = eQP(a) - eQP(i) + eh_Om(n) + dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j) + reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - ! end do ! j - ! !3p2h - ! do b=nO+1,nOrb-nR - ! num = - ERI(p,i,b,a) * eh_rho(p,b,n) * eh_rho(a,i,n) & - ! + ERI(p,i,a,b) * eh_rho(i,a,nS+n) * eh_rho(b,p,nS+n) + end do ! j + !3p2h + do b=nO+1,nOrb-nR + num = (ERI(p,i,b,a) - ERI(p,i,a,b)) * eh_rho(p,b,n) * eh_rho(a,i,n) !& + !- ERI(p,i,a,b) * eh_rho(i,a,nS+n) * eh_rho(b,p,nS+n) - ! dem1 = eQP(p) - eQP(b) - eh_Om(n) - ! dem2 = eQP(p) - eQP(b) - eQP(a) + eQP(i) - ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + dem1 = eQP(p) - eQP(b) - eh_Om(n) + dem2 = eQP(p) - eQP(b) - eQP(a) + eQP(i) + reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) & - ! - num * (reg1/dem1/dem1) * (reg2/dem2) + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) & + - num * (reg1/dem1/dem1) * (reg2/dem2) + + num = - (ERI(p,a,b,i) - ERI(p,a,i,b)) * eh_rho(p,b,n) * eh_rho(i,a,n) !& + !+ ERI(p,a,i,b) * eh_rho(a,i,nS+n) * eh_rho(b,p,nS+n) + + dem1 = eQP(a) - eQP(i) + eh_Om(n) + dem2 = eQP(p) - eQP(b) - eh_Om(n) + reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - ! ! num = - ERI(p,i,b,a) * eh_rho(p,b,n) * eh_rho(a,i,n) & - ! ! + ERI(p,i,a,b) * eh_rho(i,a,nS+n) * eh_rho(b,p,nS+n) + num = - (ERI(p,i,b,a) - ERI(p,i,a,b)) * eh_rho(p,b,nS+n) * eh_rho(a,i,nS+n) !& + !+ ERI(p,i,a,b) * eh_rho(i,a,n) * eh_rho(b,p,n) - ! ! dem1 = eQP(a) - eQP(i) - eh_Om(n) - ! ! dem2 = eQP(p) - eQP(b) - eh_Om(n) - ! ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + dem1 = eQP(a) - eQP(i) + eh_Om(n) + dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b) + reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - ! ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - ! ! num = + ERI(p,i,b,a) * eh_rho(p,b,n) * eh_rho(a,i,n) & - ! ! - ERI(p,i,a,b) * eh_rho(i,a,nS+n) * eh_rho(b,p,nS+n) - - ! ! dem1 = eQP(a) - eQP(i) - eh_Om(n) - ! ! dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b) - ! ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - - ! ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - - ! num = - ERI(p,a,b,i) * eh_rho(p,b,n) * eh_rho(i,a,n) & - ! + ERI(p,a,i,b) * eh_rho(a,i,nS+n) * eh_rho(b,p,nS+n) - - ! dem1 = eQP(a) - eQP(i) + eh_Om(n) - ! dem2 = eQP(p) - eQP(b) - eh_Om(n) - ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - - ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - - ! num = - ERI(p,i,b,a) * eh_rho(p,b,nS+n) * eh_rho(a,i,nS+n) & - ! + ERI(p,i,a,b) * eh_rho(i,a,n) * eh_rho(b,p,n) - - ! dem1 = eQP(a) - eQP(i) + eh_Om(n) - ! dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b) - ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - - ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - - ! end do ! b + end do ! b - ! end do ! n + end do ! n - ! end do ! a - ! end do ! i + end do ! a + end do ! i - ! end do ! p - ! !$OMP END DO - ! !$OMP END PARALLEL + end do ! p + !$OMP END DO + !$OMP END PARALLEL - ! call wall_time(end_t) - ! t = end_t - start_t + call wall_time(end_t) + t = end_t - start_t - ! write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building eh self-energy =',t,' seconds' - ! write(*,*) + write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building eh self-energy =',t,' seconds' + write(*,*) !-----------------------------! ! pp part of the self-energy ! @@ -290,24 +249,6 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) & - num * (reg1/dem1/dem1) * (reg2/dem2) - ! num = - ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n) - ! dem1 = hh_Om(n) - eQP(i) - eQP(j) - ! dem2 = eQP(p) + eQP(c) - hh_Om(n) - ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - - ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - - ! num = ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n) - ! dem1 = hh_Om(n) - eQP(i) - eQP(j) - ! dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j) - ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - - ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - end do ! c end do ! n end do ! j @@ -365,24 +306,6 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) & - num * (reg1/dem1/dem1) * (reg2/dem2) - ! num = ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n) - ! dem1 = ee_Om(n) - eQP(a) - eQP(b) - ! dem2 = eQP(p) + eQP(k) - ee_Om(n) - ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - - ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - - ! num = - ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n) - ! dem1 = ee_Om(n) - eQP(a) - eQP(b) - ! dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b) - ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) - - ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) - ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - end do ! c end do ! n end do ! b diff --git a/src/Parquet/G_screened_integrals.f90 b/src/Parquet/G_screened_integrals.f90 index 2e3d270..efeb6e8 100644 --- a/src/Parquet/G_screened_integrals.f90 +++ b/src/Parquet/G_screened_integrals.f90 @@ -36,15 +36,15 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb)) rho(p,q,ia) = rho(p,q,ia) & - + (ERI(q,j,p,b) - ERI(q,j,b,p)) * X !& + + (ERI(q,j,p,b) - ERI(q,j,b,p)) * X & !- (eh_Phi(q,j,b,p) + pp_Phi(q,j,p,b)) * X & - !+ (ERI(q,b,p,j) - ERI(q,b,j,p)) * Y & + + (ERI(q,b,p,j) - ERI(q,b,j,p)) * Y !& !- (eh_Phi(q,b,j,p) + pp_Phi(q,b,p,j)) * Y rho(p,q,nS+ia) = rho(p,q,nS+ia) & - + (ERI(q,b,p,j) - ERI(q,b,j,p)) * X !& + + (ERI(q,j,p,b) - ERI(q,b,j,p)) * X & !- (eh_Phi(q,j,b,p) + pp_Phi(q,j,p,b)) * X & - !+ (ERI(q,b,p,j) - ERI(q,b,j,p)) * Y & + + (ERI(q,b,p,j) - ERI(q,j,b,p)) * Y !& !- (eh_Phi(q,b,j,p) + pp_Phi(q,b,p,j)) * Y end do