From 921902962b5f482e987dd8eadb217e2e66f87359 Mon Sep 17 00:00:00 2001 From: Antoine MARIE Date: Mon, 21 Apr 2025 14:56:27 +0200 Subject: [PATCH] trying a sanity check --- src/GT/GG0T0pp.f90 | 2 +- src/GT/GGTpp_self_energy_diag.f90 | 101 ++++++- src/Parquet/G_Parquet_self_energy.f90 | 264 +++++++++++-------- src/Parquet/R_Parquet_self_energy.f90 | 362 +++++++++++++------------- 4 files changed, 434 insertions(+), 295 deletions(-) diff --git a/src/GT/GG0T0pp.f90 b/src/GT/GG0T0pp.f90 index 7be6092..2af71ca 100644 --- a/src/GT/GG0T0pp.f90 +++ b/src/GT/GG0T0pp.f90 @@ -128,7 +128,7 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T if(regularize) call GTpp_regularization(nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2) - call GGTpp_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,EcGM,Sig,Z) + call GGTpp_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,EcGM,Sig,Z,ERI) !---------------------------------------------- ! Solve the quasi-particle equation diff --git a/src/GT/GGTpp_self_energy_diag.f90 b/src/GT/GGTpp_self_energy_diag.f90 index 1f34638..f319548 100644 --- a/src/GT/GGTpp_self_energy_diag.f90 +++ b/src/GT/GGTpp_self_energy_diag.f90 @@ -1,4 +1,4 @@ -subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z) +subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z,ERI) ! Compute diagonal of the correlation part of the T-matrix self-energy @@ -20,11 +20,12 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables - integer :: i,j,a,b,p,cd,kl - double precision :: num,eps + integer :: i,j,k,a,b,c,p,m,cd,kl + double precision :: num,eps,dem1,dem2 ! Output variables @@ -72,6 +73,100 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh 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 + + ! 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 + + ! 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,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 + + ! 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 + + ! 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 + + ! 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,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 + ! end do + ! end do + ! end do + !-------------------------------------! ! Galitskii-Migdal correlation energy ! !-------------------------------------! diff --git a/src/Parquet/G_Parquet_self_energy.f90 b/src/Parquet/G_Parquet_self_energy.f90 index d47e3ff..7ab19b5 100644 --- a/src/Parquet/G_Parquet_self_energy.f90 +++ b/src/Parquet/G_Parquet_self_energy.f90 @@ -89,127 +89,151 @@ 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 - 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(j) + eh_Om(n) - reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) - reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + ! do n=1,nS + ! !3h2p + ! do j=nC+1,nO - 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(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) - 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) + ! ! 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)) - 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) - 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) + ! ! 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(j) + 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(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) + ! ! 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) - 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) * 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) * 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(b) - eh_Om(n) - 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) + ! 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,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) * 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)) + ! ! 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) + ! ! 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) + ! ! 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(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,nS+n) * eh_rho(a,i,nS+n) & - + ERI(p,i,a,b) * eh_rho(i,a,n) * eh_rho(b,p,n) + ! 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(i) - eQP(a) - eQP(b) - 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(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) + ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - end do ! b + ! 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 ! 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 ! @@ -257,22 +281,32 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& do c=nO+1,nOrb-nR 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) + dem1 = eQP(p) + eQP(c) - hh_Om(n) 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) + 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 @@ -322,22 +356,32 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& do k=nC+1,nO num = ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n) - dem1 = ee_Om(n) - eQP(a) - eQP(b) + dem1 = eQP(p) + eQP(k) - 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) + 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)) - 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) - 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 diff --git a/src/Parquet/R_Parquet_self_energy.f90 b/src/Parquet/R_Parquet_self_energy.f90 index f68a033..948771c 100644 --- a/src/Parquet/R_Parquet_self_energy.f90 +++ b/src/Parquet/R_Parquet_self_energy.f90 @@ -89,249 +89,249 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP !-------------------------------------! ! singlet eh part of the self-energy ! !-------------------------------------! - ! 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_sing_rho,eh_sing_Om,SigC,Z) -! !$OMP DO COLLAPSE(2) -! do p=nC+1,nOrb-nR + 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_sing_rho,eh_sing_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 -! num = ( - ERI(p,a,j,i) + 0.5d0*ERI(p,a,i,j))* & -! eh_sing_rho(a,i,n) * eh_sing_rho(j,p,n) + do n=1,nS + !3h2p + do j=nC+1,nO + num = ( - ERI(p,a,j,i) + 0.5d0*ERI(p,a,i,j))* & + eh_sing_rho(a,i,n) * eh_sing_rho(j,p,n) -! dem1 = eQP(a) - eQP(i) - eh_sing_Om(n) -! dem2 = eQP(p) - eQP(j) + eh_sing_Om(n) -! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) -! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + dem1 = eQP(a) - eQP(i) - eh_sing_Om(n) + dem2 = eQP(p) - eQP(j) + eh_sing_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) + 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) - 0.5d0*ERI(p,a,i,j))* & -! eh_sing_rho(a,i,n) * eh_sing_rho(j,p,n) + num = (ERI(p,a,j,i) - 0.5d0*ERI(p,a,i,j))* & + eh_sing_rho(a,i,n) * eh_sing_rho(j,p,n) -! dem1 = eQP(a) - eQP(i) - eh_sing_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_sing_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) -! num = (- ERI(p,i,j,a) + 0.5d0*ERI(p,i,a,j)) * & -! eh_sing_rho(j,p,n) * eh_sing_rho(i,a,n) + num = (- ERI(p,i,j,a) + 0.5d0*ERI(p,i,a,j)) * & + eh_sing_rho(j,p,n) * eh_sing_rho(i,a,n) -! dem1 = eQP(a) - eQP(i) + eh_sing_Om(n) -! dem2 = eQP(p) - eQP(j) + eh_sing_Om(n) -! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) -! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + dem1 = eQP(a) - eQP(i) + eh_sing_Om(n) + dem2 = eQP(p) - eQP(j) + eh_sing_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) + 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) + 0.5d0*ERI(p,a,i,j))* & -! eh_sing_rho(p,j,n) * eh_sing_rho(i,a,n) + num = (- ERI(p,a,j,i) + 0.5d0*ERI(p,a,i,j))* & + eh_sing_rho(p,j,n) * eh_sing_rho(i,a,n) -! dem1 = eQP(a) - eQP(i) + eh_sing_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_sing_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,a,b,i) + 0.5d0*ERI(p,a,i,b)) * & -! eh_sing_rho(p,b,n) * eh_sing_rho(i,a,n) + end do ! j + !3p2h + do b=nO+1,nOrb-nR + num = (- ERI(p,a,b,i) + 0.5d0*ERI(p,a,i,b)) * & + eh_sing_rho(p,b,n) * eh_sing_rho(i,a,n) -! dem1 = eQP(a) - eQP(i) + eh_sing_Om(n) -! dem2 = eQP(p) - eQP(b) - eh_sing_Om(n) -! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) -! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + dem1 = eQP(a) - eQP(i) + eh_sing_Om(n) + dem2 = eQP(p) - eQP(b) - eh_sing_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) + 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) + 0.5d0*ERI(p,i,a,b)) * & -! eh_sing_rho(b,p,n) * eh_sing_rho(i,a,n) + num = (- ERI(p,i,b,a) + 0.5d0*ERI(p,i,a,b)) * & + eh_sing_rho(b,p,n) * eh_sing_rho(i,a,n) -! dem1 = eQP(a) - eQP(i) + eh_sing_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)) + dem1 = eQP(a) - eQP(i) + eh_sing_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) + 0.5d0*ERI(p,i,a,b)) * & -! eh_sing_rho(p,b,n) * eh_sing_rho(a,i,n) + num = (- ERI(p,i,b,a) + 0.5d0*ERI(p,i,a,b)) * & + eh_sing_rho(p,b,n) * eh_sing_rho(a,i,n) -! dem1 = eQP(a) - eQP(i) - eh_sing_Om(n) -! dem2 = eQP(p) - eQP(b) - eh_sing_Om(n) -! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) -! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + dem1 = eQP(a) - eQP(i) - eh_sing_Om(n) + dem2 = eQP(p) - eQP(b) - eh_sing_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) + 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) - 0.5d0*ERI(p,i,a,b)) * & -! eh_sing_rho(p,b,n) * eh_sing_rho(a,i,n) + num = (ERI(p,i,b,a) - 0.5d0*ERI(p,i,a,b)) * & + eh_sing_rho(p,b,n) * eh_sing_rho(a,i,n) -! dem1 = eQP(a) - eQP(i) - eh_sing_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)) + dem1 = eQP(a) - eQP(i) - eh_sing_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) -! 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 -! call wall_time(end_t) -! t = end_t - start_t + end do ! p + !$OMP END DO + !$OMP END PARALLEL + call wall_time(end_t) + t = end_t - start_t -! write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for singlet eh self-energy =',t,' seconds' -! write(*,*) -! !-------------------------------------! -! ! triplet eh part of the self-energy ! -! !-------------------------------------! -! 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_trip_rho,eh_trip_Om,SigC,Z) -! !$OMP DO COLLAPSE(2) -! do p=nC+1,nOrb-nR + write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for singlet eh self-energy =',t,' seconds' + write(*,*) +!-------------------------------------! +! triplet eh part of the self-energy ! +!-------------------------------------! + 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_trip_rho,eh_trip_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 -! num = ( + 1.5d0*ERI(p,a,i,j))* & -! eh_trip_rho(a,i,n) * eh_trip_rho(j,p,n) + do n=1,nS + !3h2p + do j=nC+1,nO + num = ( + 1.5d0*ERI(p,a,i,j))* & + eh_trip_rho(a,i,n) * eh_trip_rho(j,p,n) -! dem1 = eQP(a) - eQP(i) - eh_trip_Om(n) -! dem2 = eQP(p) - eQP(j) + eh_trip_Om(n) -! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) -! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + dem1 = eQP(a) - eQP(i) - eh_trip_Om(n) + dem2 = eQP(p) - eQP(j) + eh_trip_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) + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) -! num = ( - 1.5d0*ERI(p,a,i,j))* & -! eh_trip_rho(a,i,n) * eh_trip_rho(j,p,n) + num = ( - 1.5d0*ERI(p,a,i,j))* & + eh_trip_rho(a,i,n) * eh_trip_rho(j,p,n) -! dem1 = eQP(a) - eQP(i) - eh_trip_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_trip_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) -! num = ( + 1.5d0*ERI(p,i,a,j)) * & -! eh_trip_rho(j,p,n) * eh_trip_rho(i,a,n) + num = ( + 1.5d0*ERI(p,i,a,j)) * & + eh_trip_rho(j,p,n) * eh_trip_rho(i,a,n) -! dem1 = eQP(a) - eQP(i) + eh_trip_Om(n) -! dem2 = eQP(p) - eQP(j) + eh_trip_Om(n) -! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) -! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + dem1 = eQP(a) - eQP(i) + eh_trip_Om(n) + dem2 = eQP(p) - eQP(j) + eh_trip_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) + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) -! num = ( + 1.5d0*ERI(p,a,i,j))* & -! eh_trip_rho(p,j,n) * eh_trip_rho(i,a,n) + num = ( + 1.5d0*ERI(p,a,i,j))* & + eh_trip_rho(p,j,n) * eh_trip_rho(i,a,n) -! dem1 = eQP(a) - eQP(i) + eh_trip_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_trip_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 = ( + 1.5d0*ERI(p,a,i,b)) * & -! eh_trip_rho(p,b,n) * eh_trip_rho(i,a,n) + end do ! j + !3p2h + do b=nO+1,nOrb-nR + num = ( + 1.5d0*ERI(p,a,i,b)) * & + eh_trip_rho(p,b,n) * eh_trip_rho(i,a,n) -! dem1 = eQP(a) - eQP(i) + eh_trip_Om(n) -! dem2 = eQP(p) - eQP(b) - eh_trip_Om(n) -! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) -! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + dem1 = eQP(a) - eQP(i) + eh_trip_Om(n) + dem2 = eQP(p) - eQP(b) - eh_trip_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) + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) -! num = ( + 1.5d0*ERI(p,i,a,b)) * & -! eh_trip_rho(b,p,n) * eh_trip_rho(i,a,n) + num = ( + 1.5d0*ERI(p,i,a,b)) * & + eh_trip_rho(b,p,n) * eh_trip_rho(i,a,n) -! dem1 = eQP(a) - eQP(i) + eh_trip_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)) + dem1 = eQP(a) - eQP(i) + eh_trip_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 = ( + 1.5d0*ERI(p,i,a,b)) * & -! eh_trip_rho(p,b,n) * eh_trip_rho(a,i,n) + num = ( + 1.5d0*ERI(p,i,a,b)) * & + eh_trip_rho(p,b,n) * eh_trip_rho(a,i,n) -! dem1 = eQP(a) - eQP(i) - eh_trip_Om(n) -! dem2 = eQP(p) - eQP(b) - eh_trip_Om(n) -! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) -! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) + dem1 = eQP(a) - eQP(i) - eh_trip_Om(n) + dem2 = eQP(p) - eQP(b) - eh_trip_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) + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) -! num = ( - 1.5d0*ERI(p,i,a,b)) * & -! eh_trip_rho(p,b,n) * eh_trip_rho(a,i,n) + num = ( - 1.5d0*ERI(p,i,a,b)) * & + eh_trip_rho(p,b,n) * eh_trip_rho(a,i,n) -! dem1 = eQP(a) - eQP(i) - eh_trip_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)) + dem1 = eQP(a) - eQP(i) - eh_trip_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) -! 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 -! call wall_time(end_t) -! t = end_t - start_t + end do ! p + !$OMP END DO + !$OMP END PARALLEL + call wall_time(end_t) + t = end_t - start_t -! write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for triplet eh self-energy =',t,' seconds' -! write(*,*) + write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for triplet eh self-energy =',t,' seconds' + write(*,*) !-------------------------------------! ! singlet pp part of the self-energy !