From 1d316f0fd2c4c7c93b3c4afa16e416a4a65411bd Mon Sep 17 00:00:00 2001 From: Antoine MARIE Date: Tue, 15 Apr 2025 15:55:24 +0200 Subject: [PATCH] ok with pp self-energy, eh self-energy still giving crazy results --- src/Parquet/GParquet.f90 | 4 +- src/Parquet/G_Parquet_self_energy.f90 | 212 ++++++------ src/Parquet/G_pp_Phi.f90 | 4 +- src/Parquet/RParquet.f90 | 4 +- src/Parquet/R_Parquet_self_energy.f90 | 462 +++++++++++++------------- src/Parquet/R_pp_singlet_Phi.f90 | 4 +- src/Parquet/R_pp_triplet_Phi.f90 | 4 +- 7 files changed, 348 insertions(+), 346 deletions(-) diff --git a/src/Parquet/GParquet.f90 b/src/Parquet/GParquet.f90 index 2f5e11b..5b551c0 100644 --- a/src/Parquet/GParquet.f90 +++ b/src/Parquet/GParquet.f90 @@ -8,8 +8,8 @@ subroutine GParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i ! Hard-coded parameters - logical :: print_phLR = .false. - logical :: print_ppLR = .false. + logical :: print_phLR = .true. + logical :: print_ppLR = .true. ! Input variables diff --git a/src/Parquet/G_Parquet_self_energy.f90 b/src/Parquet/G_Parquet_self_energy.f90 index 6e3f80c..e795c46 100644 --- a/src/Parquet/G_Parquet_self_energy.f90 +++ b/src/Parquet/G_Parquet_self_energy.f90 @@ -89,129 +89,129 @@ 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) + ERI(p,a,i,j))* & - eh_rho(a,i,n) * eh_rho(j,p,n) + ! do n=1,nS + ! !3h2p + ! do j=nC+1,nO + ! num = ( - ERI(p,a,j,i) + 0d0*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(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) + ! 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) - ERI(p,a,i,j))* & - eh_rho(a,i,n) * eh_rho(j,p,n) + ! num = (ERI(p,a,j,i) - 0d0*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 - 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) - num = (- ERI(p,i,j,a) + ERI(p,i,a,j)) * & - eh_rho(j,p,n) * eh_rho(i,a,n) + ! num = (- ERI(p,i,j,a) + 0d0*ERI(p,i,a,j)) * & + ! eh_rho(j,p,n) * eh_rho(i,a,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(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) + ! 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) + ERI(p,a,i,j))* & - eh_rho(p,j,n) * eh_rho(i,a,n) + ! num = (- ERI(p,a,j,i) + 0d0*ERI(p,a,i,j))* & + ! eh_rho(p,j,n) * eh_rho(i,a,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,a,b,i) + ERI(p,a,i,b)) * & - eh_rho(p,b,n) * eh_rho(i,a,n) + ! end do ! j + ! !3p2h + ! do b=nO+1,nOrb-nR + ! num = (- ERI(p,a,b,i) + 0d0*ERI(p,a,i,b)) * & + ! eh_rho(p,b,n) * eh_rho(i,a,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(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,i,b,a) + ERI(p,i,a,b)) * & - eh_rho(b,p,n) * eh_rho(i,a,n) + ! num = (- ERI(p,i,b,a) + 0d0*ERI(p,i,a,b)) * & + ! eh_rho(b,p,n) * eh_rho(i,a,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(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) + ERI(p,i,a,b)) * & - eh_rho(p,b,n) * eh_rho(a,i,n) + ! num = (- ERI(p,i,b,a) + 0d0*ERI(p,i,a,b)) * & + ! eh_rho(p,b,n) * eh_rho(a,i,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(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,i,b,a) - ERI(p,i,a,b)) * & - eh_rho(p,b,n) * eh_rho(a,i,n) + ! num = (ERI(p,i,b,a) - 0d0*ERI(p,i,a,b)) * & + ! eh_rho(p,b,n) * eh_rho(a,i,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(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 + ! 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 building eh self-energy =',t,' seconds' - write(*,*) + ! 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(*,*) + !-----------------------------! ! pp part of the self-energy ! !-----------------------------! @@ -229,7 +229,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& do n=1,nVV ! 4h1p do k=nC+1,nO - num = 2d0 * ERI(p,k,i,j) * ee_rho(i,j,n) * ee_rho(p,k,n) + num = - ERI(p,k,i,j) * ee_rho(i,j,n) * ee_rho(p,k,n) dem1 = ee_Om(n) - eQP(i) - eQP(j) dem2 = eQP(p) + eQP(k) - ee_Om(n) reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) @@ -242,7 +242,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& ! 3h2p do c=nO+1,nOrb-nR - num = 2d0 * ERI(p,c,i,j) * ee_rho(i,j,n) * ee_rho(p,c,n) + num = - ERI(p,c,i,j) * ee_rho(i,j,n) * ee_rho(p,c,n) dem1 = ee_Om(n) - eQP(i) - eQP(j) dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j) reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) @@ -257,7 +257,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& ! 3h2p do c=nO+1,nOrb-nR - num = 2d0 * ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n) + 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)) @@ -266,14 +266,14 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - num = 2d0 * ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n) + 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) + 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 @@ -294,27 +294,27 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& ! 4p1h do c=nO+1,nOrb-nR - num = 2d0 * ERI(p,c,a,b) * hh_rho(a,b,n) * hh_rho(p,c,n) + num = ERI(p,c,a,b) * hh_rho(a,b,n) * hh_rho(p,c,n) dem1 = hh_Om(n) - eQP(a) - eQP(b) 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) + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) end do ! c ! 3p2h do k=nC+1,nO - num = 2d0 * ERI(p,k,a,b) * hh_rho(a,b,n) * hh_rho(p,k,n) + num = ERI(p,k,a,b) * hh_rho(a,b,n) * hh_rho(p,k,n) dem1 = hh_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) end do ! k end do ! n @@ -322,16 +322,16 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& ! 3p2h do k=nC+1,nO - num = 2d0 * ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n) + 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) + SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) + Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - num = 2d0 * ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n) + 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)) diff --git a/src/Parquet/G_pp_Phi.f90 b/src/Parquet/G_pp_Phi.f90 index 0818b8f..1fcb813 100644 --- a/src/Parquet/G_pp_Phi.f90 +++ b/src/Parquet/G_pp_Phi.f90 @@ -31,12 +31,12 @@ subroutine G_pp_Phi(nOrb,nC,nR,nOO,nVV,ee_Om,ee_rho,hh_Om,hh_rho,pp_Phi) do n=1,nVV pp_Phi(p,q,r,s) = pp_Phi(p,q,r,s) & - + 2d0 * ee_rho(p,q,n)*ee_rho(r,s,n)/ee_Om(n) + - ee_rho(p,q,n)*ee_rho(r,s,n)/ee_Om(n) end do do n=1,nOO pp_Phi(p,q,r,s) = pp_Phi(p,q,r,s) & - - 2d0 * hh_rho(p,q,n)*hh_rho(r,s,n)/hh_Om(n) + + hh_rho(p,q,n)*hh_rho(r,s,n)/hh_Om(n) end do enddo diff --git a/src/Parquet/RParquet.f90 b/src/Parquet/RParquet.f90 index 289a89e..d7742cb 100644 --- a/src/Parquet/RParquet.f90 +++ b/src/Parquet/RParquet.f90 @@ -8,8 +8,8 @@ subroutine RParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i ! Hard-coded parameters - logical :: print_phLR = .false. - logical :: print_ppLR = .false. + logical :: print_phLR = .true. + logical :: print_ppLR = .true. ! Input variables diff --git a/src/Parquet/R_Parquet_self_energy.f90 b/src/Parquet/R_Parquet_self_energy.f90 index 45c35db..f68a033 100644 --- a/src/Parquet/R_Parquet_self_energy.f90 +++ b/src/Parquet/R_Parquet_self_energy.f90 @@ -89,249 +89,250 @@ 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 ! !-------------------------------------! @@ -347,7 +348,7 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP do n=1,nVVs ! 4h1p do k=nC+1,nO - num = ERI(p,k,i,j) * ee_sing_rho(i,j,n) * ee_sing_rho(p,k,n) + num = - 0.5d0 * ERI(p,k,i,j) * ee_sing_rho(i,j,n) * ee_sing_rho(p,k,n) dem1 = ee_sing_Om(n) - eQP(i) - eQP(j) dem2 = eQP(p) + eQP(k) - ee_sing_Om(n) reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) @@ -360,10 +361,10 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP ! 3h2p do c=nO+1,nOrb-nR - num = ERI(p,c,i,j) * ee_sing_rho(i,j,n) * ee_sing_rho(p,c,n) - !dem1 = ee_Om(n) - eQP(i) - eQP(j) + num = - 0.5d0*ERI(p,c,i,j) * ee_sing_rho(i,j,n) * ee_sing_rho(p,c,n) + dem1 = ee_sing_Om(n) - eQP(i) - eQP(j) dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j) - !reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) @@ -375,7 +376,7 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP ! 3h2p do c=nO+1,nOrb-nR - num = ERI(p,c,i,j) * hh_sing_rho(i,j,n) * hh_sing_rho(p,c,n) + num = - 0.5d0*ERI(p,c,i,j) * hh_sing_rho(i,j,n) * hh_sing_rho(p,c,n) dem1 = hh_sing_Om(n) - eQP(i) - eQP(j) dem2 = eQP(p) + eQP(c) - hh_sing_Om(n) reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) @@ -384,14 +385,14 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP 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_sing_rho(i,j,n) * hh_sing_rho(p,c,n) - !dem1 = hh_Om(n) - eQP(i) - eQP(j) + num = 0.5d0*ERI(p,c,i,j) * hh_sing_rho(i,j,n) * hh_sing_rho(p,c,n) + dem1 = hh_sing_Om(n) - eQP(i) - eQP(j) dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j) - !reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + 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 ! c end do ! n @@ -412,27 +413,27 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP ! 4p1h do c=nO+1,nOrb-nR - num = ERI(p,c,a,b) * hh_sing_rho(a,b,n) * hh_sing_rho(p,c,n) + num = 0.5d0*ERI(p,c,a,b) * hh_sing_rho(a,b,n) * hh_sing_rho(p,c,n) dem1 = hh_sing_Om(n) - eQP(a) - eQP(b) dem2 = eQP(p) + eQP(c) - hh_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) end do ! c ! 3p2h do k=nC+1,nO - num = ERI(p,k,a,b) * hh_sing_rho(a,b,n) * hh_sing_rho(p,k,n) - !dem1 = hh_Om(n) - eQP(a) - eQP(b) + num = 0.5d0*ERI(p,k,a,b) * hh_sing_rho(a,b,n) * hh_sing_rho(p,k,n) + dem1 = hh_sing_Om(n) - eQP(a) - eQP(b) dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b) - !reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + 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 ! k end do ! n @@ -440,19 +441,19 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP ! 3p2h do k=nC+1,nO - num = ERI(p,k,a,b) * ee_sing_rho(a,b,n) * ee_sing_rho(p,k,n) + num = 0.5d0*ERI(p,k,a,b) * ee_sing_rho(a,b,n) * ee_sing_rho(p,k,n) dem1 = ee_sing_Om(n) - eQP(a) - eQP(b) dem2 = eQP(p) + eQP(k) - ee_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,k,a,b) * ee_sing_rho(a,b,n) * ee_sing_rho(p,k,n) - !dem1 = ee_Om(n) - eQP(a) - eQP(b) + num = - 0.5d0*ERI(p,k,a,b) * ee_sing_rho(a,b,n) * ee_sing_rho(p,k,n) + dem1 = ee_sing_Om(n) - eQP(a) - eQP(b) dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b) - !reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) @@ -486,7 +487,7 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP do n=1,nVVt ! 4h1p do k=nC+1,nO - num = 3d0*ERI(p,k,i,j) * ee_trip_rho(i,j,n) * ee_trip_rho(p,k,n) + num = - 1.5d0 * ERI(p,k,i,j) * ee_trip_rho(i,j,n) * ee_trip_rho(p,k,n) dem1 = ee_trip_Om(n) - eQP(i) - eQP(j) dem2 = eQP(p) + eQP(k) - ee_trip_Om(n) reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) @@ -499,10 +500,10 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP ! 3h2p do c=nO+1,nOrb-nR - num = 3d0*ERI(p,c,i,j) * ee_trip_rho(i,j,n) * ee_trip_rho(p,c,n) - !dem1 = ee_Om(n) - eQP(i) - eQP(j) + num = - 1.5d0 * ERI(p,c,i,j) * ee_trip_rho(i,j,n) * ee_trip_rho(p,c,n) + dem1 = ee_trip_Om(n) - eQP(i) - eQP(j) dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j) - !reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) @@ -514,7 +515,7 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP ! 3h2p do c=nO+1,nOrb-nR - num = 3d0*ERI(p,c,i,j) * hh_trip_rho(i,j,n) * hh_trip_rho(p,c,n) + num = - 1.5d0 * ERI(p,c,i,j) * hh_trip_rho(i,j,n) * hh_trip_rho(p,c,n) dem1 = hh_trip_Om(n) - eQP(i) - eQP(j) dem2 = eQP(p) + eQP(c) - hh_trip_Om(n) reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) @@ -523,14 +524,14 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) - num = 3d0*ERI(p,c,i,j) * hh_trip_rho(i,j,n) * hh_trip_rho(p,c,n) - !dem1 = hh_Om(n) - eQP(i) - eQP(j) + num = 1.5d0 * ERI(p,c,i,j) * hh_trip_rho(i,j,n) * hh_trip_rho(p,c,n) + dem1 = hh_trip_Om(n) - eQP(i) - eQP(j) dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j) - !reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + 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 ! c end do ! n @@ -551,27 +552,27 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP ! 4p1h do c=nO+1,nOrb-nR - num = 3d0*ERI(p,c,a,b) * hh_trip_rho(a,b,n) * hh_trip_rho(p,c,n) + num = 1.5d0 * ERI(p,c,a,b) * hh_trip_rho(a,b,n) * hh_trip_rho(p,c,n) dem1 = hh_trip_Om(n) - eQP(a) - eQP(b) dem2 = eQP(p) + eQP(c) - hh_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) end do ! c ! 3p2h do k=nC+1,nO - num = 3d0*ERI(p,k,a,b) * hh_trip_rho(a,b,n) * hh_trip_rho(p,k,n) - !dem1 = hh_Om(n) - eQP(a) - eQP(b) + num = 1.5d0 * ERI(p,k,a,b) * hh_trip_rho(a,b,n) * hh_trip_rho(p,k,n) + dem1 = hh_trip_Om(n) - eQP(a) - eQP(b) dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b) - !reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + 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 ! k end do ! n @@ -579,19 +580,19 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP ! 3p2h do k=nC+1,nO - num = 3d0*ERI(p,k,a,b) * ee_trip_rho(a,b,n) * ee_trip_rho(p,k,n) + num = 1.5d0 * ERI(p,k,a,b) * ee_trip_rho(a,b,n) * ee_trip_rho(p,k,n) dem1 = ee_trip_Om(n) - eQP(a) - eQP(b) dem2 = eQP(p) + eQP(k) - ee_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 = 3d0*ERI(p,k,a,b) * ee_trip_rho(a,b,n) * ee_trip_rho(p,k,n) - !dem1 = ee_Om(n) - eQP(a) - eQP(b) + num = - 1.5d0 * ERI(p,k,a,b) * ee_trip_rho(a,b,n) * ee_trip_rho(p,k,n) + dem1 = ee_trip_Om(n) - eQP(a) - eQP(b) dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b) - !reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) + reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1)) reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2)) SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2) @@ -608,7 +609,8 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP call wall_time(end_t) t = end_t - start_t write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for triplet pp self-energy =',t,' seconds' - write(*,*) + write(*,*) + !-----------------------------! ! Renormalization factor ! !-----------------------------! diff --git a/src/Parquet/R_pp_singlet_Phi.f90 b/src/Parquet/R_pp_singlet_Phi.f90 index 73552a7..eaa275f 100644 --- a/src/Parquet/R_pp_singlet_Phi.f90 +++ b/src/Parquet/R_pp_singlet_Phi.f90 @@ -32,12 +32,12 @@ subroutine R_pp_singlet_Phi(nOrb,nC,nR,nOO,nVV,ee_sing_Om,ee_sing_rho,hh_sing_Om do n=1,nVV pp_sing_Phi(p,q,r,s) = pp_sing_Phi(p,q,r,s) & - + 2d0 * ee_sing_rho(p,q,n)*ee_sing_rho(r,s,n)/ee_sing_Om(n) + - ee_sing_rho(p,q,n)*ee_sing_rho(r,s,n)/ee_sing_Om(n) end do do n=1,nOO pp_sing_Phi(p,q,r,s) = pp_sing_Phi(p,q,r,s) & - - 2d0 * hh_sing_rho(p,q,n)*hh_sing_rho(r,s,n)/hh_sing_Om(n) + + hh_sing_rho(p,q,n)*hh_sing_rho(r,s,n)/hh_sing_Om(n) end do enddo diff --git a/src/Parquet/R_pp_triplet_Phi.f90 b/src/Parquet/R_pp_triplet_Phi.f90 index b7c204b..8d4990d 100644 --- a/src/Parquet/R_pp_triplet_Phi.f90 +++ b/src/Parquet/R_pp_triplet_Phi.f90 @@ -32,12 +32,12 @@ subroutine R_pp_triplet_Phi(nOrb,nC,nR,nOO,nVV,ee_trip_Om,ee_trip_rho,hh_trip_Om do n=1,nVV pp_trip_Phi(p,q,r,s) = pp_trip_Phi(p,q,r,s) & - + 2d0 * ee_trip_rho(p,q,n)*ee_trip_rho(r,s,n)/ee_trip_Om(n) + - ee_trip_rho(p,q,n)*ee_trip_rho(r,s,n)/ee_trip_Om(n) end do do n=1,nOO pp_trip_Phi(p,q,r,s) = pp_trip_Phi(p,q,r,s) & - - 2d0 * hh_trip_rho(p,q,n)*hh_trip_rho(r,s,n)/hh_trip_Om(n) + + hh_trip_rho(p,q,n)*hh_trip_rho(r,s,n)/hh_trip_Om(n) end do enddo