From b3859c744642d06e81c51f39221712f23390130e Mon Sep 17 00:00:00 2001 From: Antoine MARIE Date: Wed, 26 Mar 2025 16:19:23 +0100 Subject: [PATCH] saving work but not fully debugged yet --- src/Parquet/GParquet.f90 | 15 +- src/Parquet/G_eh_Phi.f90 | 37 ++++ src/Parquet/G_pp_Gam.f90 | 2 - src/Parquet/G_pp_Phi.f90 | 43 +++++ src/Parquet/G_screened_integrals.f90 | 1 + src/Parquet/RParquet.f90 | 40 +++-- src/Parquet/R_eh_singlet_Gam.f90 | 61 +------ src/Parquet/R_eh_singlet_Phi.f90 | 38 ++++ src/Parquet/R_eh_triplet_Gam.f90 | 71 +------- src/Parquet/R_eh_triplet_Phi.f90 | 38 ++++ src/Parquet/R_pp_singlet_Gam.f90 | 64 ++----- src/Parquet/R_pp_singlet_Phi.f90 | 44 +++++ src/Parquet/R_pp_triplet_Gam.f90 | 62 ++----- src/Parquet/R_pp_triplet_Phi.f90 | 44 +++++ src/Parquet/R_screened_integrals.f90 | 252 +++++++++++++++------------ 15 files changed, 457 insertions(+), 355 deletions(-) create mode 100644 src/Parquet/G_eh_Phi.f90 create mode 100644 src/Parquet/G_pp_Phi.f90 create mode 100644 src/Parquet/R_eh_singlet_Phi.f90 create mode 100644 src/Parquet/R_eh_triplet_Phi.f90 create mode 100644 src/Parquet/R_pp_singlet_Phi.f90 create mode 100644 src/Parquet/R_pp_triplet_Phi.f90 diff --git a/src/Parquet/GParquet.f90 b/src/Parquet/GParquet.f90 index 3d23baf..1fcf56d 100644 --- a/src/Parquet/GParquet.f90 +++ b/src/Parquet/GParquet.f90 @@ -9,8 +9,8 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, logical :: linearize = .true. logical :: TDA = .true. - logical :: print_phLR = .false. - logical :: print_ppLR = .false. + logical :: print_phLR = .true. + logical :: print_ppLR = .true. ! Input variables @@ -203,8 +203,9 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, end if - Aph(:,:) = Aph(:,:) + eh_Gam_A(:,:) - Bph(:,:) = Bph(:,:) + eh_Gam_B(:,:) + Aph(:,:) = Aph(:,:) + eh_Gam_A(:,:) + Bph(:,:) = Bph(:,:) + eh_Gam_B(:,:) + call phGLR(TDA,nS,Aph,Bph,EcRPA,eh_Om,XpY,XmY) @@ -258,9 +259,9 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, end if - Bpp(:,:) = Bpp(:,:) + pp_Gam_B(:,:) - Cpp(:,:) = Cpp(:,:) + pp_Gam_C(:,:) - Dpp(:,:) = Dpp(:,:) + pp_Gam_D(:,:) + ! Bpp(:,:) = Bpp(:,:) + pp_Gam_B(:,:) + ! Cpp(:,:) = Cpp(:,:) + pp_Gam_C(:,:) + ! Dpp(:,:) = Dpp(:,:) + pp_Gam_D(:,:) call ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,ee_Om,X1,Y1,hh_Om,X2,Y2,EcRPA) call wall_time(end_t) diff --git a/src/Parquet/G_eh_Phi.f90 b/src/Parquet/G_eh_Phi.f90 new file mode 100644 index 0000000..5ec86d3 --- /dev/null +++ b/src/Parquet/G_eh_Phi.f90 @@ -0,0 +1,37 @@ +subroutine G_eh_Phi(nOrb,nC,nR,nS,eh_Om,eh_rho,eh_Phi) + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nR,nS + double precision,intent(in) :: eh_Om(nS) + double precision,intent(in) :: eh_rho(nOrb,nOrb,nS+nS) + +! Local variables + integer :: p,q,r,s + integer :: n + +! Output variables + double precision, intent(out) :: eh_Phi(nOrb,nOrb,nOrb,nOrb) + +! Initialization + eh_Phi(:,:,:,:) = 0d0 + + do s = nC+1, nOrb-nR + do r = nC+1, nOrb-nR + do q = nC+1, nOrb-nR + do p = nC+1, nOrb-nR + + ! do n=1,nS + ! eh_Phi(p,q,r,s) = eh_Phi(p,q,r,s) & + ! - eh_rho(r,p,n)*eh_rho(q,s,n)/eh_Om(n) & + ! - eh_rho(p,r,nS+n)*eh_rho(s,q,nS+n)/eh_Om(n) + ! end do + + enddo + enddo + enddo + enddo + +end subroutine G_eh_Phi diff --git a/src/Parquet/G_pp_Gam.f90 b/src/Parquet/G_pp_Gam.f90 index 3b55e4a..03d99f8 100644 --- a/src/Parquet/G_pp_Gam.f90 +++ b/src/Parquet/G_pp_Gam.f90 @@ -100,8 +100,6 @@ subroutine G_pp_Gamma_B(nOrb,nC,nO,nR,nOO,nVV,eh_Phi,pp_Gam_B) ! Local variables integer :: a,b,i,j integer :: ab,ij - integer :: n - double precision,external :: Kronecker_delta ! Output variables double precision, intent(out) :: pp_Gam_B(nVV,nOO) diff --git a/src/Parquet/G_pp_Phi.f90 b/src/Parquet/G_pp_Phi.f90 new file mode 100644 index 0000000..12aad4d --- /dev/null +++ b/src/Parquet/G_pp_Phi.f90 @@ -0,0 +1,43 @@ +subroutine G_pp_Phi(nOrb,nC,nR,nOO,nVV,ee_Om,ee_rho,hh_Om,hh_rho,pp_Phi) + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nR,nOO,nVV + double precision,intent(in) :: ee_Om(nVV) + double precision,intent(in) :: ee_rho(nOrb,nOrb,nVV) + double precision,intent(in) :: hh_Om(nOO) + double precision,intent(in) :: hh_rho(nOrb,nOrb,nOO) + +! Local variables + integer :: p,q,r,s + integer :: n + +! Output variables + double precision, intent(out) :: pp_Phi(nOrb,nOrb,nOrb,nOrb) + +! Initialization + pp_Phi(:,:,:,:) = 0d0 + + do s = nC+1, nOrb-nR + do r = nC+1, nOrb-nR + do q = nC+1, nOrb-nR + do p = nC+1, nOrb-nR + + 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) + 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) + end do + + enddo + enddo + enddo + enddo + +end subroutine G_pp_Phi diff --git a/src/Parquet/G_screened_integrals.f90 b/src/Parquet/G_screened_integrals.f90 index 4d7a6b0..f6c52bf 100644 --- a/src/Parquet/G_screened_integrals.f90 +++ b/src/Parquet/G_screened_integrals.f90 @@ -146,6 +146,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 end do end do end do + end do end do ! !$OMP END DO diff --git a/src/Parquet/RParquet.f90 b/src/Parquet/RParquet.f90 index 17a79af..dbe8e42 100644 --- a/src/Parquet/RParquet.f90 +++ b/src/Parquet/RParquet.f90 @@ -141,12 +141,12 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, hh_sing_rho(:,:,:) = 0d0 hh_trip_rho(:,:,:) = 0d0 - old_eh_sing_Om(:) = 1d0 - old_eh_trip_Om(:) = 1d0 - old_ee_sing_Om(:) = 1d0 - old_ee_trip_Om(:) = 1d0 - old_hh_sing_Om(:) = 1d0 - old_hh_trip_Om(:) = 1d0 + old_eh_sing_Om(:) = 0d0 + old_eh_trip_Om(:) = 0d0 + old_ee_sing_Om(:) = 0d0 + old_ee_trip_Om(:) = 0d0 + old_hh_sing_Om(:) = 0d0 + old_hh_trip_Om(:) = 0d0 old_eh_sing_Phi(:,:,:,:) = 0d0 old_eh_trip_Phi(:,:,:,:) = 0d0 @@ -217,9 +217,11 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, old_eh_sing_Phi,old_eh_trip_Phi,old_pp_sing_Phi,old_pp_trip_Phi, & eh_sing_Gam_B) - end if + end if + Aph(:,:) = Aph(:,:) + eh_sing_Gam_A(:,:) - Bph(:,:) = Bph(:,:) + eh_sing_Gam_B(:,:) + Bph(:,:) = Bph(:,:) + eh_sing_Gam_B(:,:) + call phRLR(TDA,nS,Aph,Bph,EcRPA,eh_sing_Om,sing_XpY,sing_XmY) @@ -262,11 +264,11 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, else - call R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS, & + call R_eh_triplet_Gamma_A(nOrb,nC,nO,nR,nS, & old_eh_sing_Phi,old_eh_trip_Phi,old_pp_sing_Phi,old_pp_trip_Phi, & eh_trip_Gam_A) - if(.not.TDA) call R_eh_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS, & + if(.not.TDA) call R_eh_triplet_Gamma_B(nOrb,nC,nO,nR,nS, & old_eh_sing_Phi,old_eh_trip_Phi,old_pp_sing_Phi,old_pp_trip_Phi, & eh_trip_Gam_B) @@ -328,9 +330,9 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, end if - Bpp(:,:) = Bpp(:,:) + pp_sing_Gam_B(:,:) - Cpp(:,:) = Cpp(:,:) + pp_sing_Gam_C(:,:) - Dpp(:,:) = Dpp(:,:) + pp_sing_Gam_D(:,:) + ! Bpp(:,:) = Bpp(:,:) + pp_sing_Gam_B(:,:) + ! Cpp(:,:) = Cpp(:,:) + pp_sing_Gam_C(:,:) + ! Dpp(:,:) = Dpp(:,:) + pp_sing_Gam_D(:,:) call ppRLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,ee_sing_Om,X1s,Y1s,hh_sing_Om,X2s,Y2s,EcRPA) call wall_time(end_t) @@ -379,16 +381,16 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, else - if(.not.TDA) call R_pp_triplet_Gamma_B(nOrb,nC,nO,nR,nS,nOOt,nVVt,& + if(.not.TDA) call R_pp_triplet_Gamma_B(nOrb,nC,nO,nR,nOOt,nVVt,& old_eh_sing_Phi,old_eh_trip_Phi,pp_trip_Gam_B) - call R_pp_triplet_Gamma_C(nOrb,nO,nR,nS,nVVt,old_eh_sing_Phi,old_eh_trip_Phi,pp_trip_Gam_C) - call R_pp_triplet_Gamma_D(nOrb,nC,nO,nS,nOOt,old_eh_sing_Phi,old_eh_trip_Phi,pp_trip_Gam_D) + call R_pp_triplet_Gamma_C(nOrb,nO,nR,nVVt,old_eh_sing_Phi,old_eh_trip_Phi,pp_trip_Gam_C) + call R_pp_triplet_Gamma_D(nOrb,nC,nO,nOOt,old_eh_sing_Phi,old_eh_trip_Phi,pp_trip_Gam_D) end if - Bpp(:,:) = Bpp(:,:) + pp_trip_Gam_B(:,:) - Cpp(:,:) = Cpp(:,:) + pp_trip_Gam_C(:,:) - Dpp(:,:) = Dpp(:,:) + pp_trip_Gam_D(:,:) + ! Bpp(:,:) = Bpp(:,:) + pp_trip_Gam_B(:,:) + ! Cpp(:,:) = Cpp(:,:) + pp_trip_Gam_C(:,:) + ! Dpp(:,:) = Dpp(:,:) + pp_trip_Gam_D(:,:) call ppRLR(TDA,nOOt,nVVt,Bpp,Cpp,Dpp,ee_trip_Om,X1t,Y1t,hh_trip_Om,X2t,Y2t,EcRPA) diff --git a/src/Parquet/R_eh_singlet_Gam.f90 b/src/Parquet/R_eh_singlet_Gam.f90 index 3fe4930..7855ff6 100644 --- a/src/Parquet/R_eh_singlet_Gam.f90 +++ b/src/Parquet/R_eh_singlet_Gam.f90 @@ -14,8 +14,6 @@ subroutine R_eh_singlet_Gamma_A(nOrb,nC,nO,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing ! Local variables integer :: i,a,j,b integer :: ia,jb - integer :: n - ! Output variables double precision, intent(out) :: eh_sing_Gam_A(nS,nS) @@ -31,33 +29,8 @@ subroutine R_eh_singlet_Gamma_A(nOrb,nC,nO,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing do b=nO+1,norb-nR jb = jb + 1 - ! do n=1,nS - ! eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) & - ! + 0.5d0 * eh_sing_rho(b,a,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) & - ! + 0.5d0 * eh_sing_rho(a,b,n)*eh_sing_rho(i,j,n)/eh_sing_Om(n) & - ! + 1.5d0 * eh_trip_rho(b,a,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n) & - ! + 1.5d0 * eh_trip_rho(a,b,n)*eh_trip_rho(i,j,n)/eh_trip_Om(n) - ! end do - - ! do n=1,nVVs - ! eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) & - ! + ee_sing_rho(a,j,n)*ee_sing_rho(i,b,n)/ee_sing_Om(n) - ! end do - - ! do n=1,nOOs - ! eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) & - ! - hh_sing_rho(a,j,n)*hh_sing_rho(i,b,n)/hh_sing_Om(n) - ! end do - - ! do n=1,nVVt - ! eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) & - ! + 3d0 * ee_trip_rho(a,j,n)*ee_trip_rho(i,b,n)/ee_trip_Om(n) - ! end do - - ! do n=1,nOOt - ! eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) & - ! - 3d0 * hh_trip_rho(a,j,n)*hh_trip_rho(i,b,n)/hh_trip_Om(n) - ! end do + eh_sing_Gam_A(ia,jb) = - 0.5d0*eh_sing_Phi(a,j,b,i) - 1.5d0*eh_trip_Phi(a,j,b,i) & + + 0.5d0*pp_sing_Phi(a,j,i,b) + 1.5d0*pp_trip_Phi(a,j,i,b) enddo enddo @@ -82,7 +55,6 @@ subroutine R_eh_singlet_Gamma_B(nOrb,nC,nO,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing ! Local variables integer :: i,a,j,b integer :: ia,jb - integer :: n ! Output variables double precision, intent(out) :: eh_sing_Gam_B(nS,nS) @@ -99,33 +71,8 @@ subroutine R_eh_singlet_Gamma_B(nOrb,nC,nO,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing do b=nO+1,norb-nR jb = jb + 1 - ! do n=1,nS - ! eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) & - ! + 0.5d0 * eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) & - ! + 0.5d0 * eh_sing_rho(a,j,n)*eh_sing_rho(i,b,n)/eh_sing_Om(n) & - ! + 1.5d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) & - ! + 1.5d0 * eh_trip_rho(a,j,n)*eh_trip_rho(i,b,n)/eh_trip_Om(n) - ! end do - - ! do n=1,nVVs - ! eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) & - ! + ee_sing_rho(a,b,n)*ee_sing_rho(i,j,n)/ee_sing_Om(n) - ! end do - - ! do n=1,nOOs - ! eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) & - ! - hh_sing_rho(a,b,n)*hh_sing_rho(i,j,n)/hh_sing_Om(n) - ! end do - - ! do n=1,nVVt - ! eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) & - ! + 3d0 * ee_trip_rho(a,b,n)*ee_trip_rho(i,j,n)/ee_trip_Om(n) - ! end do - - ! do n=1,nOOt - ! eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) & - ! - 3d0 * hh_trip_rho(a,b,n)*hh_trip_rho(i,j,n)/hh_trip_Om(n) - ! end do + eh_sing_Gam_B(ia,jb) = - 0.5d0*eh_sing_Phi(a,b,j,i) - 1.5d0*eh_trip_Phi(a,b,j,i) & + + 0.5d0*pp_sing_Phi(a,b,i,j) + 1.5d0*pp_trip_Phi(a,b,i,j) enddo enddo diff --git a/src/Parquet/R_eh_singlet_Phi.f90 b/src/Parquet/R_eh_singlet_Phi.f90 new file mode 100644 index 0000000..5971e26 --- /dev/null +++ b/src/Parquet/R_eh_singlet_Phi.f90 @@ -0,0 +1,38 @@ +subroutine R_eh_singlet_Phi(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_sing_Phi) + + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nR,nS + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS+nS) + +! Local variables + integer :: p,q,r,s + integer :: n + +! Output variables + double precision,intent(out) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb) + +! Initialization + eh_sing_Phi(:,:,:,:) = 0d0 + + do s = nC+1, nOrb-nR + do r = nC+1, nOrb-nR + do q = nC+1, nOrb-nR + do p = nC+1, nOrb-nR + + ! do n=1,nS + ! eh_sing_Phi(p,q,r,s) = eh_sing_Phi(p,q,r,s) & + ! - eh_sing_rho(r,p,n)*eh_sing_rho(q,s,n)/eh_sing_Om(n) & + ! - eh_sing_rho(p,r,n+nS)*eh_sing_rho(s,q,n+nS)/eh_sing_Om(n) + ! end do + + enddo + enddo + enddo + enddo + +end subroutine R_eh_singlet_Phi diff --git a/src/Parquet/R_eh_triplet_Gam.f90 b/src/Parquet/R_eh_triplet_Gam.f90 index 796b810..a0e9ae3 100644 --- a/src/Parquet/R_eh_triplet_Gam.f90 +++ b/src/Parquet/R_eh_triplet_Gam.f90 @@ -1,11 +1,11 @@ -subroutine R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing_Phi,pp_trip_Phi,eh_trip_Gam_A) +subroutine R_eh_triplet_Gamma_A(nOrb,nC,nO,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing_Phi,pp_trip_Phi,eh_trip_Gam_A) ! Compute irreducible vertex in the triplet pp channel implicit none ! Input variables - integer,intent(in) :: nOrb,nC,nO,nV,nR,nS + integer,intent(in) :: nOrb,nC,nO,nR,nS double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: pp_sing_Phi(nOrb,nOrb,nOrb,nOrb) @@ -14,8 +14,6 @@ subroutine R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_s ! Local variables integer :: i,a,j,b integer :: ia,jb - integer :: n - double precision,external :: Kronecker_delta ! Output variables double precision, intent(out) :: eh_trip_Gam_A(nS,nS) @@ -32,33 +30,8 @@ subroutine R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_s do b=nO+1,norb-nR jb = jb + 1 - ! do n=1,nS - ! eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) & - ! + 0.5d0 * eh_sing_rho(b,a,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) & - ! + 0.5d0 * eh_sing_rho(a,b,n)*eh_sing_rho(i,j,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_trip_rho(b,a,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n) & - ! - 0.5d0 * eh_trip_rho(a,b,n)*eh_trip_rho(i,j,n)/eh_trip_Om(n) - ! end do - - ! do n=1,nVVs - ! eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) & - ! - ee_sing_rho(a,j,n)*ee_sing_rho(i,b,n)/ee_sing_Om(n) - ! end do - - ! do n=1,nOOs - ! eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) & - ! + hh_sing_rho(a,j,n)*hh_sing_rho(i,b,n)/hh_sing_Om(n) - ! end do - - ! do n=1,nVVt - ! eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) & - ! + ee_trip_rho(a,j,n)*ee_trip_rho(i,b,n)/ee_trip_Om(n) - ! end do - - ! do n=1,nOOt - ! eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) & - ! - hh_trip_rho(a,j,n)*hh_trip_rho(i,b,n)/hh_trip_Om(n) - ! end do + eh_trip_Gam_A(ia,jb) = - 0.5d0*eh_sing_Phi(a,j,b,i) + 0.5d0*eh_trip_Phi(a,j,b,i) & + - 0.5d0*pp_sing_Phi(a,j,i,b) + 0.5d0*pp_trip_Phi(a,j,i,b) enddo enddo @@ -67,14 +40,13 @@ subroutine R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_s end subroutine R_eh_triplet_Gamma_A -subroutine R_eh_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing_Phi,pp_trip_Phi,eh_trip_Gam_B) - +subroutine R_eh_triplet_Gamma_B(nOrb,nC,nO,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing_Phi,pp_trip_Phi,eh_trip_Gam_B) ! Compute irreducible vertex in the triplet pp channel implicit none ! Input variables - integer,intent(in) :: nOrb,nC,nO,nV,nR,nS + integer,intent(in) :: nOrb,nC,nO,nR,nS double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: pp_sing_Phi(nOrb,nOrb,nOrb,nOrb) @@ -83,8 +55,6 @@ subroutine R_eh_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_s ! Local variables integer :: i,a,j,b integer :: ia,jb - integer :: n - double precision,external :: Kronecker_delta ! Output variables double precision, intent(out) :: eh_trip_Gam_B(nS,nS) @@ -101,33 +71,8 @@ subroutine R_eh_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_s do b=nO+1,norb-nR jb = jb + 1 - ! do n=1,nS - ! eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) & - ! + 0.5d0 * eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) & - ! + 0.5d0 * eh_sing_rho(a,j,n)*eh_sing_rho(i,b,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) & - ! - 0.5d0 * eh_trip_rho(a,j,n)*eh_trip_rho(i,b,n)/eh_trip_Om(n) - ! end do - - ! do n=1,nVVs - ! eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) & - ! - ee_sing_rho(a,b,n)*ee_sing_rho(i,j,n)/ee_sing_Om(n) - ! end do - - ! do n=1,nOOs - ! eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) & - ! + hh_sing_rho(a,b,n)*hh_sing_rho(i,j,n)/hh_sing_Om(n) - ! end do - - ! do n=1,nVVt - ! eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) & - ! + ee_trip_rho(a,b,n)*ee_trip_rho(i,j,n)/ee_trip_Om(n) - ! end do - - ! do n=1,nOOt - ! eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) & - ! - hh_trip_rho(a,b,n)*hh_trip_rho(i,j,n)/hh_trip_Om(n) - ! end do + eh_trip_Gam_B(ia,jb) = - 0.5d0*eh_sing_Phi(a,b,j,i) + 0.5d0*eh_trip_Phi(a,b,j,i) & + - 0.5d0*pp_sing_Phi(a,b,i,j) + 0.5d0*pp_trip_Phi(a,b,i,j) enddo enddo diff --git a/src/Parquet/R_eh_triplet_Phi.f90 b/src/Parquet/R_eh_triplet_Phi.f90 new file mode 100644 index 0000000..a6ec43f --- /dev/null +++ b/src/Parquet/R_eh_triplet_Phi.f90 @@ -0,0 +1,38 @@ +subroutine R_eh_triplet_Phi(nOrb,nC,nR,nS,eh_trip_Om,eh_trip_rho,eh_trip_Phi) + + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nR,nS + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS+nS) + +! Local variables + integer :: p,q,r,s + integer :: n + +! Output variables + double precision,intent(out) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb) + +! Initialization + eh_trip_Phi(:,:,:,:) = 0d0 + + do s = nC+1, nOrb-nR + do r = nC+1, nOrb-nR + do q = nC+1, nOrb-nR + do p = nC+1, nOrb-nR + + ! do n=1,nS + ! eh_trip_Phi(p,q,r,s) = eh_trip_Phi(p,q,r,s) & + ! - eh_trip_rho(r,p,n)*eh_trip_rho(q,s,n)/eh_trip_Om(n) & + ! - eh_trip_rho(p,r,n+nS)*eh_trip_rho(s,q,n+nS)/eh_trip_Om(n) + ! end do + + enddo + enddo + enddo + enddo + +end subroutine R_eh_triplet_Phi diff --git a/src/Parquet/R_pp_singlet_Gam.f90 b/src/Parquet/R_pp_singlet_Gam.f90 index 2a9fe2f..cd9987c 100644 --- a/src/Parquet/R_pp_singlet_Gam.f90 +++ b/src/Parquet/R_pp_singlet_Gam.f90 @@ -1,17 +1,17 @@ -subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nS,nOOs,eh_sing_Phi,eh_trip_Phi,pp_sing_Gam_D) +subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nOOs,eh_sing_Phi,eh_trip_Phi,pp_sing_Gam_D) ! Compute irreducible vertex in the triplet pp channel implicit none ! Input variables - integer,intent(in) :: nOrb,nC,nO,nS,nOOs + integer,intent(in) :: nOrb,nC,nO,nOOs double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb) ! Local variables integer :: i,j,k,l integer :: ij,kl - integer :: n + double precision,external :: Kronecker_delta ! Output variables double precision, intent(out) :: pp_sing_Gam_D(nOOs,nOOs) @@ -34,19 +34,10 @@ subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nS,nOOs,eh_sing_Phi,eh_trip_Phi,pp_si do l=k,nO kl = kl +1 - ! do n=1,nS - ! pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl) & - ! - 0.5d0 * eh_sing_rho(k,i,n)*eh_sing_rho(j,l,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_sing_rho(i,k,n)*eh_sing_rho(l,j,n)/eh_sing_Om(n) & - ! + 1.5d0 * eh_trip_rho(k,i,n)*eh_trip_rho(j,l,n)/eh_trip_Om(n) & - ! + 1.5d0 * eh_trip_rho(i,k,n)*eh_trip_rho(l,j,n)/eh_trip_Om(n) & - ! - 0.5d0 * eh_sing_rho(l,i,n)*eh_sing_rho(j,k,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_sing_rho(i,l,n)*eh_sing_rho(k,j,n)/eh_sing_Om(n) & - ! + 1.5d0 * eh_trip_rho(l,i,n)*eh_trip_rho(j,k,n)/eh_trip_Om(n) & - ! + 1.5d0 * eh_trip_rho(i,l,n)*eh_trip_rho(k,j,n)/eh_trip_Om(n) - ! end do + pp_sing_Gam_D(ij,kl) = 0.5d0*eh_sing_Phi(i,j,k,l) - 1.5d0*eh_trip_Phi(i,j,k,l) & + - 1.5d0*eh_sing_Phi(i,j,l,k) + 0.5d0*eh_trip_Phi(i,j,l,k) - ! pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) end do end do @@ -57,20 +48,20 @@ subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nS,nOOs,eh_sing_Phi,eh_trip_Phi,pp_si end subroutine -subroutine R_pp_singlet_Gamma_C(nOrb,nO,nR,nS,nVVs,eh_sing_Phi,eh_trip_Phi,pp_sing_Gam_C) +subroutine R_pp_singlet_Gamma_C(nOrb,nO,nR,nVVs,eh_sing_Phi,eh_trip_Phi,pp_sing_Gam_C) ! Compute irreducible vertex in the triplet pp channel implicit none ! Input variables - integer,intent(in) :: nOrb,nO,nR,nS,nVVs + integer,intent(in) :: nOrb,nO,nR,nVVs double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb) ! Local variables integer :: a,b,c,d integer :: ab,cd - integer :: n + double precision,external :: Kronecker_delta ! Output variables double precision, intent(out) :: pp_sing_Gam_C(nVVs,nVVs) @@ -93,19 +84,11 @@ subroutine R_pp_singlet_Gamma_C(nOrb,nO,nR,nS,nVVs,eh_sing_Phi,eh_trip_Phi,pp_si do d=c,nOrb - nR cd = cd +1 - ! do n=1,nS - ! pp_sing_Gam_C(ab,cd) = pp_sing_Gam_C(ab,cd) & - ! - 0.5d0 * eh_sing_rho(c,a,n)*eh_sing_rho(b,d,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_sing_rho(a,c,n)*eh_sing_rho(d,b,n)/eh_sing_Om(n) & - ! + 1.5d0 * eh_trip_rho(c,a,n)*eh_trip_rho(b,d,n)/eh_trip_Om(n) & - ! + 1.5d0 * eh_trip_rho(a,c,n)*eh_trip_rho(d,b,n)/eh_trip_Om(n) & - ! - 0.5d0 * eh_sing_rho(d,a,n)*eh_sing_rho(b,c,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_sing_rho(a,d,n)*eh_sing_rho(c,b,n)/eh_sing_Om(n) & - ! + 1.5d0 * eh_trip_rho(d,a,n)*eh_trip_rho(b,c,n)/eh_trip_Om(n) & - ! + 1.5d0 * eh_trip_rho(a,d,n)*eh_trip_rho(c,b,n)/eh_trip_Om(n) - ! end do - ! pp_sing_Gam_C(ab,cd) = pp_sing_Gam_C(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + pp_sing_Gam_C(ab,cd) = 0.5d0*eh_sing_Phi(a,b,c,d) - 1.5d0*eh_trip_Phi(a,b,c,d) & + - 1.5d0*eh_sing_Phi(a,b,d,c) + 0.5d0*eh_trip_Phi(a,b,d,c) + + pp_sing_Gam_C(ab,cd) = pp_sing_Gam_C(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) end do end do @@ -116,20 +99,20 @@ subroutine R_pp_singlet_Gamma_C(nOrb,nO,nR,nS,nVVs,eh_sing_Phi,eh_trip_Phi,pp_si end subroutine -subroutine R_pp_singlet_Gamma_B(nOrb,nC,nO,nR,nS,nOOs,nVVs,eh_sing_Phi,eh_trip_Phi,pp_sing_Gam_B) +subroutine R_pp_singlet_Gamma_B(nOrb,nC,nO,nR,nOOs,nVVs,eh_sing_Phi,eh_trip_Phi,pp_sing_Gam_B) ! Compute irreducible vertex in the triplet pp channel implicit none ! Input variables - integer,intent(in) :: nOrb,nC,nO,nR,nS,nOOs,nVVs + integer,intent(in) :: nOrb,nC,nO,nR,nOOs,nVVs double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb) ! Local variables integer :: a,b,i,j integer :: ab,ij - integer :: n + double precision,external :: Kronecker_delta ! Output variables double precision, intent(out) :: pp_sing_Gam_B(nVVs,nOOs) @@ -152,19 +135,10 @@ subroutine R_pp_singlet_Gamma_B(nOrb,nC,nO,nR,nS,nOOs,nVVs,eh_sing_Phi,eh_trip_P do j=i,nO ij = ij +1 - ! do n=1,nS - ! pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij) & - ! - 0.5d0 * eh_sing_rho(i,a,n)*eh_sing_rho(b,j,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_sing_rho(a,i,n)*eh_sing_rho(j,b,n)/eh_sing_Om(n) & - ! + 1.5d0 * eh_trip_rho(i,a,n)*eh_trip_rho(b,j,n)/eh_trip_Om(n) & - ! + 1.5d0 * eh_trip_rho(a,i,n)*eh_trip_rho(j,b,n)/eh_trip_Om(n) & - ! - 0.5d0 * eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_sing_rho(a,j,n)*eh_sing_rho(i,b,n)/eh_sing_Om(n) & - ! + 1.5d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) & - ! + 1.5d0 * eh_trip_rho(a,j,n)*eh_trip_rho(i,b,n)/eh_trip_Om(n) - ! end do + pp_sing_Gam_B(ab,ij) = 0.5d0*eh_sing_Phi(a,b,i,j) - 1.5d0*eh_trip_Phi(a,b,i,j) & + - 1.5d0*eh_sing_Phi(a,b,j,i) + 0.5d0*eh_trip_Phi(a,b,j,i) - ! pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) + pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) end do end do diff --git a/src/Parquet/R_pp_singlet_Phi.f90 b/src/Parquet/R_pp_singlet_Phi.f90 new file mode 100644 index 0000000..9a0f245 --- /dev/null +++ b/src/Parquet/R_pp_singlet_Phi.f90 @@ -0,0 +1,44 @@ +subroutine R_pp_singlet_Phi(nOrb,nC,nR,nOO,nVV,ee_sing_Om,ee_sing_rho,hh_sing_Om,hh_sing_rho,pp_sing_Phi) + + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nR,nOO,nVV + double precision,intent(in) :: ee_sing_Om(nVV) + double precision,intent(in) :: ee_sing_rho(nOrb,nOrb,nVV) + double precision,intent(in) :: hh_sing_Om(nOO) + double precision,intent(in) :: hh_sing_rho(nOrb,nOrb,nOO) + +! Local variables + integer :: p,q,r,s + integer :: n + +! Output variables + double precision,intent(out) :: pp_sing_Phi(nOrb,nOrb,nOrb,nOrb) + +! Initialization + pp_sing_Phi(:,:,:,:) = 0d0 + + do s = nC+1, nOrb-nR + do r = nC+1, nOrb-nR + do q = nC+1, nOrb-nR + do p = nC+1, nOrb-nR + + do n=1,nVV + pp_sing_Phi(p,q,r,s) = pp_sing_Phi(p,q,r,s) & + + 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) & + - hh_sing_rho(p,q,n)*hh_sing_rho(r,s,n)/hh_sing_Om(n) + end do + + enddo + enddo + enddo + enddo + +end subroutine R_pp_singlet_Phi diff --git a/src/Parquet/R_pp_triplet_Gam.f90 b/src/Parquet/R_pp_triplet_Gam.f90 index 973589f..75bdc24 100644 --- a/src/Parquet/R_pp_triplet_Gam.f90 +++ b/src/Parquet/R_pp_triplet_Gam.f90 @@ -1,10 +1,10 @@ -subroutine R_pp_triplet_Gamma_D(nOrb,nC,nO,nS,nOOt,eh_sing_Phi,eh_trip_Phi,pp_trip_Gam_D) +subroutine R_pp_triplet_Gamma_D(nOrb,nC,nO,nOOt,eh_sing_Phi,eh_trip_Phi,pp_trip_Gam_D) ! Compute irreducible vertex in the triplet pp channel implicit none ! Input variables - integer,intent(in) :: nOrb,nC,nO,nS,nOOt + integer,intent(in) :: nOrb,nC,nO,nOOt double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb) @@ -34,20 +34,9 @@ subroutine R_pp_triplet_Gamma_D(nOrb,nC,nO,nS,nOOt,eh_sing_Phi,eh_trip_Phi,pp_tr do l=k+1,nO kl = kl +1 - do n=1,nS - - ! pp_trip_Gam_D(ij,kl) = pp_trip_Gam_D(ij,kl) & - ! - 0.5d0 * eh_sing_rho(k,i,n)*eh_sing_rho(j,l,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_sing_rho(i,k,n)*eh_sing_rho(l,j,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_trip_rho(k,i,n)*eh_trip_rho(j,l,n)/eh_trip_Om(n) & - ! - 0.5d0 * eh_trip_rho(i,k,n)*eh_trip_rho(l,j,n)/eh_trip_Om(n) & - ! + 0.5d0 * eh_sing_rho(l,i,n)*eh_sing_rho(j,k,n)/eh_sing_Om(n) & - ! + 0.5d0 * eh_sing_rho(i,l,n)*eh_sing_rho(k,j,n)/eh_sing_Om(n) & - ! + 0.5d0 * eh_trip_rho(l,i,n)*eh_trip_rho(j,k,n)/eh_trip_Om(n) & - ! + 0.5d0 * eh_trip_rho(i,l,n)*eh_trip_rho(k,j,n)/eh_trip_Om(n) - - end do - + pp_trip_Gam_D(ij,kl) = 0.5d0*eh_sing_Phi(i,j,k,l) + 0.5d0*eh_trip_Phi(i,j,k,l) & + - 0.5d0*eh_sing_Phi(i,j,l,k) - 0.5d0*eh_trip_Phi(i,j,l,k) + end do end do end do @@ -57,13 +46,13 @@ subroutine R_pp_triplet_Gamma_D(nOrb,nC,nO,nS,nOOt,eh_sing_Phi,eh_trip_Phi,pp_tr end subroutine -subroutine R_pp_triplet_Gamma_C(nOrb,nO,nR,nS,nVVt,eh_sing_Phi,eh_trip_Phi,pp_trip_Gam_C) +subroutine R_pp_triplet_Gamma_C(nOrb,nO,nR,nVVt,eh_sing_Phi,eh_trip_Phi,pp_trip_Gam_C) ! Compute irreducible vertex in the triplet pp channel implicit none ! Input variables - integer,intent(in) :: nOrb,nO,nR,nS,nVVt + integer,intent(in) :: nOrb,nO,nR,nVVt double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb) @@ -92,23 +81,13 @@ subroutine R_pp_triplet_Gamma_C(nOrb,nO,nR,nS,nVVt,eh_sing_Phi,eh_trip_Phi,pp_tr do c=nO+1,nOrb - nR do d=c+1,nOrb - nR cd = cd +1 - - do n=1,nS - - ! pp_trip_Gam_C(ab,cd) = pp_trip_Gam_C(ab,cd) & - ! - 0.5d0 * eh_sing_rho(c,a,n)*eh_sing_rho(b,d,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_sing_rho(a,c,n)*eh_sing_rho(d,b,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_trip_rho(c,a,n)*eh_trip_rho(b,d,n)/eh_trip_Om(n) & - ! - 0.5d0 * eh_trip_rho(a,c,n)*eh_trip_rho(d,b,n)/eh_trip_Om(n) & - ! + 0.5d0 * eh_sing_rho(d,a,n)*eh_sing_rho(b,c,n)/eh_sing_Om(n) & - ! + 0.5d0 * eh_sing_rho(a,d,n)*eh_sing_rho(c,b,n)/eh_sing_Om(n) & - ! + 0.5d0 * eh_trip_rho(d,a,n)*eh_trip_rho(b,c,n)/eh_trip_Om(n) & - ! + 0.5d0 * eh_trip_rho(a,d,n)*eh_trip_rho(c,b,n)/eh_trip_Om(n) - - end do + + pp_trip_Gam_C(ab,cd) = 0.5d0*eh_sing_Phi(a,b,c,d) + 0.5d0*eh_trip_Phi(a,b,c,d) & + - 0.5d0*eh_sing_Phi(a,b,d,c) - 0.5d0*eh_trip_Phi(a,b,d,c) end do end do + end do end do ! !$OMP END DO @@ -116,13 +95,13 @@ subroutine R_pp_triplet_Gamma_C(nOrb,nO,nR,nS,nVVt,eh_sing_Phi,eh_trip_Phi,pp_tr end subroutine -subroutine R_pp_triplet_Gamma_B(nOrb,nC,nO,nR,nS,nOOt,nVVt,eh_sing_Phi,eh_trip_Phi,pp_trip_Gam_B) +subroutine R_pp_triplet_Gamma_B(nOrb,nC,nO,nR,nOOt,nVVt,eh_sing_Phi,eh_trip_Phi,pp_trip_Gam_B) ! Compute irreducible vertex in the triplet pp channel implicit none ! Input variables - integer,intent(in) :: nOrb,nC,nO,nR,nS,nOOt,nVVt + integer,intent(in) :: nOrb,nC,nO,nR,nOOt,nVVt double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb) @@ -152,19 +131,8 @@ subroutine R_pp_triplet_Gamma_B(nOrb,nC,nO,nR,nS,nOOt,nVVt,eh_sing_Phi,eh_trip_P do j=i+1,nO ij = ij +1 - do n=1,nS - - ! pp_trip_Gam_B(ab,ij) = pp_trip_Gam_B(ab,ij) & - ! - 0.5d0 * eh_sing_rho(i,a,n)*eh_sing_rho(b,j,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_sing_rho(a,i,n)*eh_sing_rho(j,b,n)/eh_sing_Om(n) & - ! - 0.5d0 * eh_trip_rho(i,a,n)*eh_trip_rho(b,j,n)/eh_trip_Om(n) & - ! - 0.5d0 * eh_trip_rho(a,i,n)*eh_trip_rho(j,b,n)/eh_trip_Om(n) & - ! + 0.5d0 * eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) & - ! + 0.5d0 * eh_sing_rho(a,j,n)*eh_sing_rho(i,b,n)/eh_sing_Om(n) & - ! + 0.5d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) & - ! + 0.5d0 * eh_trip_rho(a,j,n)*eh_trip_rho(i,b,n)/eh_trip_Om(n) - - end do + pp_trip_Gam_B(ab,ij) = 0.5d0*eh_sing_Phi(a,b,i,j) + 0.5d0*eh_trip_Phi(a,b,i,j) & + - 0.5d0*eh_sing_Phi(a,b,j,i) - 0.5d0*eh_trip_Phi(a,b,j,i) end do end do diff --git a/src/Parquet/R_pp_triplet_Phi.f90 b/src/Parquet/R_pp_triplet_Phi.f90 new file mode 100644 index 0000000..5af9613 --- /dev/null +++ b/src/Parquet/R_pp_triplet_Phi.f90 @@ -0,0 +1,44 @@ +subroutine R_pp_triplet_Phi(nOrb,nC,nR,nOO,nVV,ee_trip_Om,ee_trip_rho,hh_trip_Om,hh_trip_rho,pp_trip_Phi) + + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nR,nOO,nVV + double precision,intent(in) :: ee_trip_Om(nVV) + double precision,intent(in) :: ee_trip_rho(nOrb,nOrb,nVV) + double precision,intent(in) :: hh_trip_Om(nOO) + double precision,intent(in) :: hh_trip_rho(nOrb,nOrb,nOO) + +! Local variables + integer :: p,q,r,s + integer :: n + +! Output variables + double precision,intent(out) :: pp_trip_Phi(nOrb,nOrb,nOrb,nOrb) + +! Initialization + pp_trip_Phi(:,:,:,:) = 0d0 + + do s = nC+1, nOrb-nR + do r = nC+1, nOrb-nR + do q = nC+1, nOrb-nR + do p = nC+1, nOrb-nR + + do n=1,nVV + pp_trip_Phi(p,q,r,s) = pp_trip_Phi(p,q,r,s) & + + 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) & + - hh_trip_rho(p,q,n)*hh_trip_rho(r,s,n)/hh_trip_Om(n) + end do + + enddo + enddo + enddo + enddo + +end subroutine R_pp_triplet_Phi diff --git a/src/Parquet/R_screened_integrals.f90 b/src/Parquet/R_screened_integrals.f90 index c9c5a74..4e2ad64 100644 --- a/src/Parquet/R_screened_integrals.f90 +++ b/src/Parquet/R_screened_integrals.f90 @@ -17,7 +17,7 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr double precision :: X,Y ! Output variables - double precision,intent(out) :: rho(nOrb,nOrb,nS) + double precision,intent(out) :: rho(nOrb,nOrb,nS+nS) rho(:,:,:) = 0d0 ! !$OMP PARALLEL & @@ -33,18 +33,20 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr do b=nO+1,nOrb-nR jb = jb + 1 - ! do ia=1,nS + do ia=1,nS - ! X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb)) - ! Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb)) + X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb)) + Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb)) - ! rho(p,q,ia) = rho(p,q,ia) & - ! + (2d0*ERI(p,j,q,b) - ERI(p,j,b,q))*X & - ! + (2d0*ERI(p,b,q,j) - ERI(p,b,j,q))*Y & - ! + 1d0*eh_sing_Gam(p,j,q,b)*X & - ! + 1d0*eh_sing_Gam(p,b,q,j)*Y + rho(p,q,ia) = (2d0*ERI(q,j,p,b) - ERI(q,j,b,p) & + - 0.5d0*eh_sing_Phi(q,j,b,p) - 1.5d0*eh_trip_Phi(q,j,b,p) & + + 0.5d0*pp_sing_Phi(q,j,p,b) + 1.5d0*pp_trip_Phi(q,j,p,b)) * X - ! end do + rho(p,q,nS+ia) = (2d0*ERI(q,b,p,j) - ERI(q,b,j,p) & + - 0.5d0*eh_sing_Phi(q,b,j,p) - 1.5d0*eh_trip_Phi(q,b,j,p) & + + 0.5d0*pp_sing_Phi(q,b,p,j) + 1.5d0*pp_trip_Phi(q,b,p,j)) * Y + + end do end do end do @@ -75,7 +77,7 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr double precision :: X,Y ! Output variables - double precision,intent(out) :: rho(nOrb,nOrb,nS) + double precision,intent(out) :: rho(nOrb,nOrb,nS+nS) rho(:,:,:) = 0d0 ! !$OMP PARALLEL & @@ -85,24 +87,30 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr ! !$OMP DO do q=nC+1,nOrb-nR do p=nC+1,nOrb-nR + jb = 0 do j=nC+1,nO do b=nO+1,nOrb-nR jb = jb + 1 + do ia=1,nS - ! X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb)) - ! Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb)) + X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb)) + Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb)) - ! rho(p,q,ia) = rho(p,q,ia) & - ! - ERI(p,j,b,q)*X & - ! - ERI(p,b,j,q)*Y & - ! + 1d0*eh_trip_Gam(p,j,q,b)*X & - ! + 1d0*eh_trip_Gam(p,b,q,j)*Y + rho(p,q,ia) = (- ERI(q,j,b,p) & + - 0.5d0*eh_sing_Phi(q,j,b,p) + 0.5d0*eh_trip_Phi(q,j,b,p) & + - 0.5d0*pp_sing_Phi(q,j,p,b) + 0.5d0*pp_trip_Phi(q,j,p,b)) * X + + rho(p,q,nS+ia) = (- ERI(q,b,j,p) & + - 0.5d0*eh_sing_Phi(q,b,j,p) + 0.5d0*eh_trip_Phi(q,b,j,p) & + - 0.5d0*pp_sing_Phi(q,b,p,j) + 0.5d0*pp_trip_Phi(q,b,p,j)) * Y end do + end do end do + end do end do ! !$OMP END DO @@ -156,62 +164,71 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi, ! !$OMP DO COLLAPSE(2) do q=nC+1,nOrb-nR do p=nC+1,nOrb-nR - - ! ab = 0 - ! do a=nO+1,nOrb-nR - ! do b=a,nOrb-nR - ! ab = ab + 1 - ! cd = 0 - ! do c=nO+1,nOrb-nR - ! do d=c,nOrb-nR - ! cd = cd + 1 - ! rho1(p,q,ab) = rho1(p,q,ab) & - ! + (ERI(p,q,c,d) + ERI(p,q,d,c) + 1d0*pp_sing_Gam(p,q,c,d))*X1(cd,ab) & - ! /sqrt(1d0 + Kronecker_delta(c,d)) - ! end do - ! end do - - ! kl = 0 - ! do k=nC+1,nO - ! do l=k,nO - ! kl = kl + 1 - ! rho1(p,q,ab) = rho1(p,q,ab) & - ! + (ERI(p,q,k,l) + ERI(p,q,l,k) + 1d0*pp_sing_Gam(p,q,k,l))*Y1(kl,ab) & - ! /sqrt(1d0 + Kronecker_delta(k,l)) - ! end do - ! end do + ab=0 + do a = nO+1, nOrb-nR + do b = a, nOrb-nR + ab = ab + 1 - ! end do - ! end do - - ! ij = 0 - ! do i=nC+1,nO - ! do j=i,nO - ! ij = ij + 1 + cd = 0 + do c = nO+1, nOrb-nR + do d = c, nOrb-nR + cd = cd + 1 + + rho1(p,q,ab) = (ERI(p,q,c,d) + ERI(p,q,d,c) & + + 0.5d0*eh_sing_Phi(p,q,c,d) - 1.5d0*eh_trip_Phi(p,q,c,d) & + - 1.5d0*eh_sing_Phi(p,q,d,c) + 0.5d0*eh_trip_Phi(p,q,d,c))& + *X1(cd,ab)/sqrt(1d0 + Kronecker_delta(c,d)) + + end do ! d + end do ! c + + kl = 0 + do k = nC+1, nO + do l = k, nO + + kl = kl + 1 + + rho1(p,q,ab) = (ERI(p,q,k,l) + ERI(p,q,l,k) & + + 0.5d0*eh_sing_Phi(p,q,k,l) - 1.5d0*eh_trip_Phi(p,q,k,l) & + - 1.5d0*eh_sing_Phi(p,q,l,k) + 0.5d0*eh_trip_Phi(p,q,l,k))& + *Y1(kl,ab)/sqrt(1d0 + Kronecker_delta(k,l)) + end do ! l + end do ! k + end do ! b + end do ! a + + ij = 0 + do i = nC+1, nO + do j = i, nO + ij = ij + 1 + + cd = 0 + do c = nO+1, nOrb-nR + do d = c, nOrb-nR + cd = cd + 1 + + rho2(p,q,ij) = (ERI(p,q,c,d) + ERI(p,q,d,c) & + + 0.5d0*eh_sing_Phi(p,q,c,d) - 1.5d0*eh_trip_Phi(p,q,c,d) & + - 1.5d0*eh_sing_Phi(p,q,d,c) + 0.5d0*eh_trip_Phi(p,q,d,c))& + *X2(cd,ij)/sqrt(1d0 + Kronecker_delta(c,d)) + end do ! d + end do ! c + + kl = 0 + do k = nC+1, nO + do l = k, nO + kl = kl + 1 + + rho2(p,q,ij) = (ERI(p,q,k,l) + ERI(p,q,l,k) & + + 0.5d0*eh_sing_Phi(p,q,k,l) - 1.5d0*eh_trip_Phi(p,q,k,l) & + - 1.5d0*eh_sing_Phi(p,q,l,k) + 0.5d0*eh_trip_Phi(p,q,l,k))& + *Y2(kl,ij)/sqrt(1d0 + Kronecker_delta(k,l)) + end do ! l + end do ! k - ! cd = 0 - ! do c=nO+1,nOrb-nR - ! do d=c,nOrb-nR - ! cd = cd + 1 - ! rho2(p,q,ij) = rho2(p,q,ij) & - ! + (ERI(p,q,c,d) + ERI(p,q,d,c) + 1d0*pp_sing_Gam(p,q,c,d))*X2(cd,ij) & - ! /sqrt(1d0 + Kronecker_delta(c,d)) - ! end do - ! end do - - ! kl = 0 - ! do k=nC+1,nO - ! do l=k,nO - ! kl = kl + 1 - ! rho2(p,q,ij) = rho2(p,q,ij) & - ! + (ERI(p,q,k,l) + ERI(p,q,l,k) + 1d0*pp_sing_Gam(p,q,k,l))*Y2(kl,ij) & - ! /sqrt(1d0 + Kronecker_delta(k,l)) - ! end do - ! end do - - ! end do - ! end do + end do ! j + end do ! i end do end do @@ -264,62 +281,67 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi, ! !$OMP DO COLLAPSE(2) do q = nC+1, nOrb-nR do p = nC+1, nOrb-nR + ab = 0 - - ! do a = nO+1, nOrb-nR - ! do b = a+1, nOrb-nR - ! ab = ab + 1 + do a = nO+1, nOrb-nR + do b = a+1, nOrb-nR + ab = ab + 1 - ! cd = 0 - ! do c = nO+1, nOrb-nR - ! do d = c+1, nOrb-nR - ! cd = cd + 1 + cd = 0 + do c = nO+1, nOrb-nR + do d = c+1, nOrb-nR + cd = cd + 1 - ! rho1(p,q,ab) = rho1(p,q,ab) & - ! + (ERI(p,q,c,d) - ERI(p,q,d,c) + 1d0*pp_trip_Gam(p,q,c,d))*X1(cd,ab) - ! end do ! d - ! end do ! c + rho1(p,q,ab) = (ERI(p,q,c,d) - ERI(p,q,d,c) & + + 0.5d0*eh_sing_Phi(p,q,c,d) + 0.5d0*eh_trip_Phi(p,q,c,d) & + - 0.5d0*eh_sing_Phi(p,q,d,c) - 0.5d0*eh_trip_Phi(p,q,d,c) )*X1(cd,ab) + + end do ! d + end do ! c - ! kl = 0 - ! do k = nC+1, nO - ! do l = k+1, nO + kl = 0 + do k = nC+1, nO + do l = k+1, nO - ! kl = kl + 1 + kl = kl + 1 - ! rho1(p,q,ab) = rho1(p,q,ab) & - ! + (ERI(p,q,k,l) - ERI(p,q,l,k) + 1d0*pp_trip_Gam(p,q,k,l))*Y1(kl,ab) - ! end do ! l - ! end do ! k - ! end do ! b - ! end do ! a + rho1(p,q,ab) = (ERI(p,q,k,l) - ERI(p,q,l,k) & + + 0.5d0*eh_sing_Phi(p,q,k,l) + 0.5d0*eh_trip_Phi(p,q,k,l) & + - 0.5d0*eh_sing_Phi(p,q,l,k) - 0.5d0*eh_trip_Phi(p,q,l,k) )*Y1(kl,ab) + end do ! l + end do ! k + end do ! b + end do ! a - ! ij = 0 - ! do i = nC+1, nO - ! do j = i+1, nO - ! ij = ij + 1 + ij = 0 + do i = nC+1, nO + do j = i+1, nO + ij = ij + 1 - ! cd = 0 - ! do c = nO+1, nOrb-nR - ! do d = c+1, nOrb-nR - ! cd = cd + 1 + cd = 0 + do c = nO+1, nOrb-nR + do d = c+1, nOrb-nR + cd = cd + 1 - ! rho2(p,q,ij) = rho2(p,q,ij) & - ! + (ERI(p,q,c,d) - ERI(p,q,d,c) + 1d0*pp_trip_Gam(p,q,c,d))*X2(cd,ij) - ! end do ! d - ! end do ! c + rho2(p,q,ij) = (ERI(p,q,c,d) - ERI(p,q,d,c) & + + 0.5d0*eh_sing_Phi(p,q,c,d) + 0.5d0*eh_trip_Phi(p,q,c,d) & + - 0.5d0*eh_sing_Phi(p,q,d,c) - 0.5d0*eh_trip_Phi(p,q,d,c) )*X2(cd,ij) + end do ! d + end do ! c - ! kl = 0 - ! do k = nC+1, nO - ! do l = k+1, nO - ! kl = kl + 1 + kl = 0 + do k = nC+1, nO + do l = k+1, nO + kl = kl + 1 - ! rho2(p,q,ij) = rho2(p,q,ij) & - ! + (ERI(p,q,k,l) - ERI(p,q,l,k) + 1d0*pp_trip_Gam(p,q,k,l))*Y2(kl,ij) - ! end do ! l - ! end do ! k + rho2(p,q,ij) = (ERI(p,q,k,l) - ERI(p,q,l,k) & + + 0.5d0*eh_sing_Phi(p,q,k,l) + 0.5d0*eh_trip_Phi(p,q,k,l) & + - 0.5d0*eh_sing_Phi(p,q,l,k) - 0.5d0*eh_trip_Phi(p,q,l,k) )*Y2(kl,ij) + end do ! l + end do ! k - ! end do ! j - ! end do ! i + end do ! j + end do ! i end do ! p end do ! q