diff --git a/src/Parquet/GParquet.f90 b/src/Parquet/GParquet.f90 index d8d26b9..c4e534c 100644 --- a/src/Parquet/GParquet.f90 +++ b/src/Parquet/GParquet.f90 @@ -497,8 +497,8 @@ subroutine GParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i eOld(:) = eQP(:) ! Print for one-body part - - call print_parquet_1b(nOrb,nO,eHF,SigC,eQP,Z,n_it_1b,err_1b,ENuc,EGHF,EcGM,Ec_eh,Ec_pp) + + call G_print_parquet_1b(nOrb,nO,eHF,SigC,eQP,Z,n_it_1b,err_1b,ENuc,EGHF,EcGM,Ec_eh,Ec_pp) deallocate(eQPlin,Z,SigC) diff --git a/src/Parquet/G_Parquet_self_energy.f90 b/src/Parquet/G_Parquet_self_energy.f90 new file mode 100644 index 0000000..783e5cf --- /dev/null +++ b/src/Parquet/G_Parquet_self_energy.f90 @@ -0,0 +1,357 @@ +subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& + eh_rho,eh_Om,ee_rho,ee_Om,hh_rho,hh_Om,EcGM,SigC,Z) + +! Compute correlation part of the self-energy coming from irreducible vertices contribution + implicit none + include 'parameters.h' + +! Input variables + double precision,intent(in) :: eta + integer,intent(in) :: nOrb + integer,intent(in) :: nC, nO, nV, nR + integer,intent(in) :: nS, nOO, nVV + double precision,intent(in) :: eQP(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eh_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_Om(nS) + double precision,intent(in) :: ee_rho(nOrb,nOrb,nVV) + double precision,intent(in) :: ee_Om(nVV) + double precision,intent(in) :: hh_rho(nOrb,nOrb,nOO) + double precision,intent(in) :: hh_Om(nOO) + +! Local variables + integer :: i,j,k,a,b,c + integer :: p,n + double precision :: eps,dem1,dem2,reg,reg1,reg2 + double precision :: num + double precision :: start_t,end_t,t + +! Output variables + double precision,intent(out) :: SigC(nOrb) + double precision,intent(out) :: Z(nOrb) + double precision,intent(out) :: EcGM + + ! Initialize + + SigC(:) = 0d0 + Z(:) = 0d0 + EcGM = 0d0 + +!-----------------------------! +! GF2 part of the self-energy ! +!-----------------------------! + call wall_time(start_t) + do p=nC+1,nOrb-nR + ! 2h1p sum + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nOrb-nR + + 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 + + SigC(p) = SigC(p) + num*reg/eps + Z(p) = Z(p) - num*reg/eps**2 + + end do + end do + end do + ! 2p1h sum + do i=nC+1,nO + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR + + 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 + + SigC(p) = SigC(p) + num*reg/eps + Z(p) = Z(p) - num*reg/eps**2 + + end do + end do + end do + end do + call wall_time(end_t) + t = end_t - start_t + + write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building GF(2) self-energy =',t,' seconds' + write(*,*) +! !-----------------------------! +! ! 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_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 n=1,nS +! !3h2p +! do j=nC+1,nO +! num = ERI(p,a,j,i) * & +! (eh_rho(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,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(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,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)) + +! 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(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,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)) + +! 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(j,p,n) * eh_rho(a,i,n) - eh_rho(j,i,n) * eh_rho(a,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) + + +! end do ! j +! !3p2h +! do b=nO+1,nOrb-nR +! num = ERI(p,a,b,i) * & +! (eh_rho(b,p,n) * eh_rho(i,a,n) - eh_rho(b,a,n) * eh_rho(i,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)) + +! 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(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,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) + +! !num = ERI(p,i,b,a) * & +! !(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,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) + +! !num = ERI(p,i,b,a) * & +! !(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,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)) + +! 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 ! a +! end do ! i + +! 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(*,*) +! !-----------------------------! +! ! pp part of the self-energy ! +! !-----------------------------! +! call wall_time(start_t) +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(p,i,j,k,c,n,num,dem1,dem2,reg1,reg2) & +! !$OMP SHARED(nC,nO,nOrb,nR,nOO,nVV,eta,ERI,eQP,ee_rho,ee_Om,hh_rho,hh_Om,SigC,Z) +! !$OMP DO COLLAPSE(2) +! do p=nC+1,nOrb-nR + +! do i=nC+1,nO +! do j=nC+1,nO +! do n=1,nVV +! ! 4h1p +! do k=nC+1,nO +! 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)) +! 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 ! k +! ! 3h2p +! do c=nO+1,nOrb-nR + +! 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)) +! 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 ! a +! end do ! n +! do n=1,nOO +! ! 3h2p +! 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) +! 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 +! end do ! i + +! end do ! p +! !$OMP END DO +! !$OMP END PARALLEL +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(p,k,a,b,c,n,num,dem1,dem2,reg1,reg2) & +! !$OMP SHARED(nC,nO,nOrb,nR,nOO,nVV,eta,ERI,eQP,ee_rho,ee_Om,hh_rho,hh_Om,SigC,Z) +! !$OMP DO COLLAPSE(2) +! do p=nC+1,nOrb-nR +! do a=nO+1,nOrb-nR +! do b=nO+1,nOrb-nR +! do n=1,nOO +! ! 4p1h +! do c=nO+1,nOrb-nR + +! 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) + +! end do ! c +! ! 3p2h +! do k=nC+1,nO + +! 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) + +! end do ! k +! end do ! n +! do n=1,nVV +! ! 3p2h +! 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) +! 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 +! end do ! a + +! 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 pp self-energy =',t,' seconds' +! write(*,*) + +!-----------------------------! +! Renormalization factor ! +!-----------------------------! + + call vecout(nOrb,Z) + Z(:) = 1d0/(1d0 - Z(:)) + +!-------------------------------------! +! Galitskii-Migdal correlation energy ! +!-------------------------------------! + + EcGM = 0d0 + +end subroutine diff --git a/src/Parquet/G_irred_Parquet_self_energy.f90 b/src/Parquet/G_irred_Parquet_self_energy.f90 deleted file mode 100644 index 2a3b143..0000000 --- a/src/Parquet/G_irred_Parquet_self_energy.f90 +++ /dev/null @@ -1,349 +0,0 @@ -subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& - eh_rho,eh_Om,ee_rho,ee_Om,hh_rho,hh_Om,EcGM,SigC,Z) - -! Compute correlation part of the self-energy coming from irreducible vertices contribution - implicit none - include 'parameters.h' - -! Input variables - double precision,intent(in) :: eta - integer,intent(in) :: nOrb - integer,intent(in) :: nC, nO, nV, nR - integer,intent(in) :: nS, nOO, nVV - double precision,intent(in) :: eQP(nOrb) - double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) - double precision,intent(in) :: eh_rho(nOrb,nOrb,nS) - double precision,intent(in) :: eh_Om(nS) - double precision,intent(in) :: ee_rho(nOrb,nOrb,nVV) - double precision,intent(in) :: ee_Om(nVV) - double precision,intent(in) :: hh_rho(nOrb,nOrb,nOO) - double precision,intent(in) :: hh_Om(nOO) - -! Local variables - integer :: i,j,k,a,b,c - integer :: p,n - double precision :: eps,dem1,dem2,reg,reg1,reg2 - double precision :: num - double precision :: start_t,end_t,t - -! Output variables - double precision,intent(out) :: SigC(nOrb) - double precision,intent(out) :: Z(nOrb) - double precision,intent(out) :: EcGM - - ! Initialize - - SigC(:) = 0d0 - Z(:) = 0d0 - EcGM = 0d0 - -!-----------------------------! -! GF2 part of the self-energy ! -!-----------------------------! - call wall_time(start_t) - do p=nC+1,nOrb-nR - ! 2h1p sum - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nOrb-nR - - 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 - - SigC(p) = SigC(p) + num*reg/eps - Z(p) = Z(p) - num*reg/eps**2 - - end do - end do - end do - ! 2p1h sum - do i=nC+1,nO - do a=nO+1,nOrb-nR - do b=nO+1,nOrb-nR - - 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 - - SigC(p) = SigC(p) + num*reg/eps - Z(p) = Z(p) - num*reg/eps**2 - - end do - end do - end do - end do - call wall_time(end_t) - t = end_t - start_t - - write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building GF(2) self-energy =',t,' seconds' - write(*,*) -!-----------------------------! -! 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_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 n=1,nS - !3h2p - do j=nC+1,nO - num = ERI(p,a,j,i) * & - (eh_rho(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,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(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,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)) - - 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(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,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)) - - 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(j,p,n) * eh_rho(a,i,n) - eh_rho(j,i,n) * eh_rho(a,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) - - - end do ! j - !3p2h - do b=nO+1,nOrb-nR - num = ERI(p,a,b,i) * & - (eh_rho(b,p,n) * eh_rho(i,a,n) - eh_rho(b,a,n) * eh_rho(i,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)) - - 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(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,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) - - !num = ERI(p,i,b,a) * & - !(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,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) - - !num = ERI(p,i,b,a) * & - !(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,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)) - - 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 ! a - end do ! i - - 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(*,*) -!-----------------------------! -! pp part of the self-energy ! -!-----------------------------! - call wall_time(start_t) - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(p,i,j,k,c,n,num,dem1,dem2,reg1,reg2) & - !$OMP SHARED(nC,nO,nOrb,nR,nOO,nVV,eta,ERI,eQP,ee_rho,ee_Om,hh_rho,hh_Om,SigC,Z) - !$OMP DO COLLAPSE(2) - do p=nC+1,nOrb-nR - - do i=nC+1,nO - do j=nC+1,nO - do n=1,nVV - ! 4h1p - do k=nC+1,nO - 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)) - 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 ! k - ! 3h2p - do c=nO+1,nOrb-nR - - 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)) - 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 ! a - end do ! n - do n=1,nOO - ! 3h2p - 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) - 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 - end do ! i - - end do ! p - !$OMP END DO - !$OMP END PARALLEL - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(p,k,a,b,c,n,num,dem1,dem2,reg1,reg2) & - !$OMP SHARED(nC,nO,nOrb,nR,nOO,nVV,eta,ERI,eQP,ee_rho,ee_Om,hh_rho,hh_Om,SigC,Z) - !$OMP DO COLLAPSE(2) - do p=nC+1,nOrb-nR - do a=nO+1,nOrb-nR - do b=nO+1,nOrb-nR - do n=1,nOO - ! 4p1h - do c=nO+1,nOrb-nR - - 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) - - end do ! c - ! 3p2h - do k=nC+1,nO - - 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) - - end do ! k - end do ! n - do n=1,nVV - ! 3p2h - 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) - 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 - end do ! a - - 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 pp self-energy =',t,' seconds' - write(*,*) -!-----------------------------! -! Renormalization factor ! -!-----------------------------! - - Z(:) = 1d0/(1d0 - Z(:)) - -end subroutine diff --git a/src/Parquet/RParquet.f90 b/src/Parquet/RParquet.f90 index 310e4cc..a228c7c 100644 --- a/src/Parquet/RParquet.f90 +++ b/src/Parquet/RParquet.f90 @@ -600,7 +600,7 @@ subroutine RParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i write(*,*)' Two-body convergence failed ' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - stop + !stop else @@ -624,13 +624,17 @@ subroutine RParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i write(*,*) 'Building self-energy' call wall_time(start_t) - call R_irred_Parquet_self_energy(nOrb,nC,nO,nV,nR,eOld,EcGM,SigC,Z) + call R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eOld,ERI, & + eh_sing_rho,old_eh_sing_Om,eh_trip_rho,old_eh_trip_Om, & + ee_sing_rho,old_ee_sing_Om,ee_trip_rho,old_ee_trip_Om, & + hh_sing_rho,old_hh_sing_Om,hh_trip_rho,old_hh_trip_Om, & + EcGM,SigC,Z) call wall_time(end_t) t = end_t - start_t write(*,'(A50,1X,F9.3,A8)') 'Wall time for self energy =',t,' seconds' write(*,*) - eQPlin(:) = eHF(:) !+ Z(:)*SigC(:) + eQPlin(:) = eHF(:) + Z(:)*SigC(:) ! Solve the quasi-particle equation @@ -647,15 +651,19 @@ subroutine RParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i write(*,*) stop - end if - - deallocate(eQPlin,Z,SigC) - + end if + ! Check one-body converge err_1b = maxval(abs(eOld - eQP)) eOld(:) = eQP(:) - + + ! Print for one-body part + + call R_print_parquet_1b(nOrb,nO,eHF,SigC,eQP,Z,n_it_1b,err_1b,ENuc,ERHF,EcGM,Ec_eh,Ec_pp) + + deallocate(eQPlin,Z,SigC) + call wall_time(end_1b) t_1b = end_1b - start_1b write(*,'(A50,1X,F9.3,A8)') 'Wall time for one-body iteration =',t_1b,' seconds' diff --git a/src/Parquet/R_Parquet_self_energy.f90 b/src/Parquet/R_Parquet_self_energy.f90 new file mode 100644 index 0000000..1e85c31 --- /dev/null +++ b/src/Parquet/R_Parquet_self_energy.f90 @@ -0,0 +1,102 @@ +subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP,ERI, & + eh_sing_rho,eh_sing_Om,eh_trip_rho,eh_trip_Om, & + ee_sing_rho,ee_sing_Om,ee_trip_rho,ee_trip_Om, & + hh_sing_rho,hh_sing_Om,hh_trip_rho,hh_trip_Om, & + EcGM,SigC,Z) + +! Compute correlation part of the self-energy with only irreducible vertices contribution + implicit none + include 'parameters.h' + +! Input variables + double precision,intent(in) :: eta + integer,intent(in) :: nOrb,nC,nO,nV,nR + integer,intent(in) :: nS,nOOs,nVVs,nOOt,nVVt + double precision,intent(in) :: eQP(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: ee_sing_rho(nOrb,nOrb,nVVs) + double precision,intent(in) :: ee_sing_Om(nVVs) + double precision,intent(in) :: ee_trip_rho(nOrb,nOrb,nVVt) + double precision,intent(in) :: ee_trip_Om(nVVt) + double precision,intent(in) :: hh_sing_rho(nOrb,nOrb,nOOs) + double precision,intent(in) :: hh_sing_Om(nOOs) + double precision,intent(in) :: hh_trip_rho(nOrb,nOrb,nOOt) + double precision,intent(in) :: hh_trip_Om(nOOt) + +! Local variables + integer :: i,j,k,a,b,c + integer :: p,n + double precision :: eps,dem1,dem2,reg,reg1,reg2 + double precision :: num + double precision :: start_t,end_t,t + +! Output variables + double precision,intent(out) :: EcGM + double precision,intent(out) :: SigC(nOrb) + double precision,intent(out) :: Z(nOrb) + +! Initialize + + SigC(:) = 0d0 + Z(:) = 0d0 + EcGM = 0d0 + +!-----------------------------! +! GF2 part of the self-energy ! +!-----------------------------! + call wall_time(start_t) + do p=nC+1,nOrb-nR + ! 2h1p sum + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nOrb-nR + + eps = eQP(p) + eQP(a) - eQP(i) - eQP(j) + reg = (1d0 - exp(- 2d0 * eta * eps * eps)) + num = ERI(p,a,j,i)*(2d0*ERI(j,i,p,a) - ERI(j,i,a,p)) + + SigC(p) = SigC(p) + num*reg/eps + Z(p) = Z(p) - num*reg/eps**2 + + end do + end do + end do + ! 2p1h sum + do i=nC+1,nO + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR + + eps = eQP(p) + eQP(i) - eQP(a) - eQP(b) + reg = (1d0 - exp(- 2d0 * eta * eps * eps)) + num = ERI(p,i,b,a)*(2d0*ERI(b,a,p,i) - ERI(b,a,i,p)) + + SigC(p) = SigC(p) + num*reg/eps + Z(p) = Z(p) - num*reg/eps**2 + + end do + end do + end do + end do + call wall_time(end_t) + t = end_t - start_t + + write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building GF(2) self-energy =',t,' seconds' + write(*,*) + +!-----------------------------! +! Renormalization factor ! +!-----------------------------! + call vecout(nOrb,Z) + Z(:) = 1d0/(1d0 - Z(:)) + +!-------------------------------------! +! Galitskii-Migdal correlation energy ! +!-------------------------------------! + + EcGM = 0d0 + +end subroutine diff --git a/src/Parquet/R_irred_Parquet_self_energy.f90 b/src/Parquet/R_irred_Parquet_self_energy.f90 deleted file mode 100644 index d0a0bbd..0000000 --- a/src/Parquet/R_irred_Parquet_self_energy.f90 +++ /dev/null @@ -1,94 +0,0 @@ -subroutine R_irred_Parquet_self_energy(nOrb,nC,nO,nV,nR,e,EcGM,SigC,Z) - -! Compute correlation part of the self-energy with only irreducible vertices contribution - implicit none - include 'parameters.h' - -! Input variables - integer,intent(in) :: nOrb,nC,nO,nV,nR - double precision,intent(in) :: e(nOrb) -! Local variables - integer :: p,i,j,a,b - double precision :: D2p1h,D2h1p -! Output variables - double precision,intent(out) :: EcGM - double precision,intent(out) :: SigC(nOrb) - double precision,intent(out) :: Z(nOrb) - -!----------------------------! -! Static Parquet self-energy ! -!----------------------------! - SigC(:) = 0d0 - ! 2h1p part of the correlation self-energy - do p=nC+1,nOrb-nR - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nOrb-nR - - D2h1p = e(p) + e(a) - e(i) - e(j) - SigC(p) = SigC(p) !+ 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim - - end do - end do - end do - end do - ! 2p1h part of the correlation self-energy - do p=nC+1,nOrb-nR - do i=nC+1,nO - do a=nO+1,nOrb-nR - do b=nO+1,nOrb-nR - - D2p1h = e(p) + e(i) - e(a) - e(b) - SigC(p) = SigC(p) !+ 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam - - end do - end do - end do - end do -!------------------------! -! Renormalization factor ! -!------------------------! - Z(:) = 0d0 - ! 2h1p part of the renormlization factor - do p=nC+1,nOrb-nR - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nOrb-nR - - D2h1p = e(p) + e(a) - e(i) - e(j) - Z(p) = Z(p) - - end do - end do - end do - end do - ! 2p1h part of the renormlization factor - do p=nC+1,nOrb-nR - do i=nC+1,nO - do a=nO+1,nOrb-nR - do b=nO+1,nOrb-nR - - D2p1h = e(p) + e(i) - e(a) - e(b) - Z(p) = Z(p) - - end do - end do - end do - end do - - Z(:) = 1d0/(1d0 - Z(:)) - -!-------------------------------------! -! Galitskii-Migdal correlation energy ! -!-------------------------------------! - - EcGM = 0d0 - ! do i=nC+1,nO - ! do a=nO+1,nOrb-nR - ! do m=1,nS - - ! end do - ! end do - ! end do - -end subroutine diff --git a/src/Parquet/print_parquet_1b.f90 b/src/Parquet/print_parquet_1b.f90 deleted file mode 100644 index a2324d9..0000000 --- a/src/Parquet/print_parquet_1b.f90 +++ /dev/null @@ -1,64 +0,0 @@ -subroutine print_parquet_1b(nOrb,nO,eHF,SigC,eQP,Z,n_it_1b,err_1b,ENuc,EGHF,EcGM,Ec_eh,Ec_pp) - -! Print one-electron energies and other stuff for G0F2 - - implicit none - include 'parameters.h' - - integer,intent(in) :: nOrb - integer,intent(in) :: nO - double precision,intent(in) :: eHF(nOrb) - double precision,intent(in) :: SigC(nOrb) - double precision,intent(in) :: eQP(nOrb) - double precision,intent(in) :: Z(nOrb) - integer,intent(in) :: n_it_1b - double precision,intent(in) :: err_1b - double precision,intent(in) :: ENuc - double precision,intent(in) :: EGHF - double precision,intent(in) :: EcGM - double precision,intent(in) :: Ec_eh - double precision,intent(in) :: Ec_pp - - integer :: p - integer :: HOMO - integer :: LUMO - double precision :: Gap - -! HOMO and LUMO - - HOMO = nO - LUMO = HOMO + 1 - Gap = eQP(LUMO) - eQP(HOMO) - -! Dump results - - write(*,*)'-------------------------------------------------------------------------------' - write(*,*)' Parquet self-energy ' - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & - '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' - write(*,*)'-------------------------------------------------------------------------------' - - do p=1,nOrb - write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & - '|',p,'|',eHF(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eQP(p)*HaToeV,'|' - end do - - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A60,I15)') 'One-body iteration # ',n_it_1b - write(*,'(2X,A60,F15.6)') 'One-body convergence ',err_1b - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A60,F15.6,A3)') 'Parquet HOMO energy = ',eQP(HOMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'Parquet LUMO energy = ',eQP(LUMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'Parquet HOMO-LUMO gap = ',Gap*HaToeV,' eV' - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A60,F15.6,A3)') ' Parquet total energy = ',ENuc + EGHF + EcGM,' au' - write(*,'(2X,A60,F15.6,A3)') ' Parquet correlation energy = ',EcGM,' au' - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A60,F15.6,A3)') ' eh-RPA correlation energy = ',Ec_eh,' au' - write(*,'(2X,A60,F15.6,A3)') ' pp-RPA correlation energy = ',Ec_pp,' au' - write(*,'(2X,A60,F15.6,A3)') '(eh+pp)-RPA correlation energy = ',Ec_pp,' au' - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) - -end subroutine