10
1
mirror of https://github.com/pfloos/quack synced 2025-04-01 06:21:37 +02:00

saving work but still not debug

This commit is contained in:
Antoine Marie 2025-03-26 22:47:31 +01:00
parent ae2fd5c9ab
commit 9a571c5b10
8 changed files with 109 additions and 110 deletions

View File

@ -117,7 +117,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
! Memory allocation
allocate(old_eh_Om(nS),old_ee_Om(nVV),old_hh_Om(nOO))
allocate(eh_rho(nOrb,nOrb,nS+nS),ee_rho(nOrb,nOrb,nVV),hh_rho(nOrb,nOrb,nOO))
allocate(eh_rho(nOrb,nOrb,nS),ee_rho(nOrb,nOrb,nVV),hh_rho(nOrb,nOrb,nOO))
allocate(old_eh_Phi(nOrb,nOrb,nOrb,nOrb),old_pp_Phi(nOrb,nOrb,nOrb,nOrb))
! Initialization
@ -259,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)
@ -328,7 +328,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
deallocate(eh_rho,ee_rho,hh_rho)
! TODO Once we will compute the blocks of kernel starting from the 4-tensors we can move the freeing up
! Memory allocation
allocate(eh_rho(nOrb,nOrb,nS+nS))
allocate(eh_rho(nOrb,nOrb,nS))
allocate(ee_rho(nOrb,nOrb,nVV),hh_rho(nOrb,nOrb,nOO))
! Build singlet eh integrals
@ -384,7 +384,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for pp reducible kernel =',t,' seconds'
write(*,*)
! alpha = 0.01d0
! alpha = 0.1d0
! eh_Phi(:,:,:,:) = alpha * eh_Phi(:,:,:,:) + (1d0 - alpha) * old_eh_Phi(:,:,:,:)
! pp_Phi(:,:,:,:) = alpha * pp_Phi(:,:,:,:) + (1d0 - alpha) * old_pp_Phi(:,:,:,:)

View File

@ -1,12 +1,12 @@
subroutine G_eh_Phi(nOrb,nC,nR,nS,eh_Om,eh_rho,eh_Phi)
! Compute irreducible vertex in the triplet pp channel
! Compute irreducible vertex in the eh 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)
double precision,intent(in) :: eh_rho(nOrb,nOrb,nS)
! Local variables
integer :: p,q,r,s
@ -23,11 +23,11 @@ subroutine G_eh_Phi(nOrb,nC,nR,nS,eh_Om,eh_rho,eh_Phi)
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
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,n)*eh_rho(s,q,n)/eh_Om(n)
end do
enddo
enddo

View File

@ -15,7 +15,7 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho
double precision :: X,Y
! Output variables
double precision,intent(out) :: rho(nOrb,nOrb,nS+nS)
double precision,intent(out) :: rho(nOrb,nOrb,nS)
rho(:,:,:) = 0d0
! !$OMP PARALLEL &
@ -25,26 +25,28 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho
! !$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
do ia=1,nS
jb = 0
do j=nC+1,nO
do b=nO+1,nOrb-nR
jb = jb + 1
X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb))
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
rho(p,q,ia) = (ERI(q,j,p,b) - ERI(q,j,b,p)) * X &
+ (- eh_Phi(q,j,b,p) + pp_Phi(q,j,p,b)) * X
rho(p,q,nS+ia) = (ERI(q,b,p,j) - ERI(q,b,j,p)) * Y &
+ (- eh_Phi(q,b,j,p) + pp_Phi(q,b,p,j)) * Y
rho(p,q,ia) = rho(p,q,ia) &
+ (ERI(q,j,p,b) - ERI(q,j,b,p) &
- 0d0*eh_Phi(q,j,b,p) + 0d0*pp_Phi(q,j,p,b)) * X &
+ (ERI(q,b,p,j) - ERI(q,b,j,p) &
- 0d0*eh_Phi(q,b,j,p) + 0d0*pp_Phi(q,b,p,j)) * Y
end do
end do
end do
end do
end do
! !$OMP END DO
@ -76,7 +78,6 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
integer :: a,b,c,d
integer :: p,q
integer :: ab,cd,ij,kl
double precision,external :: Kronecker_delta
! Output variables
@ -106,8 +107,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
do c=nO+1,nOrb-nR
do d=c+1,nOrb-nR
cd = cd + 1
rho1(p,q,ab) = ( ERI(p,q,c,d) - ERI(p,q,d,c) ) * X1(cd,ab) &
+ ( eh_Phi(p,q,c,d) - eh_Phi(p,q,d,c) ) * X1(cd,ab)
rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,c,d) - ERI(p,q,d,c) &
+ 0d0*eh_Phi(p,q,c,d) - 0d0*eh_Phi(p,q,d,c) ) * X1(cd,ab)
end do
end do
@ -115,8 +116,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho1(p,q,ab) = ( ERI(p,q,k,l) - ERI(p,q,l,k) ) * Y1(kl,ab) &
+ ( eh_Phi(p,q,k,l) - eh_Phi(p,q,l,k) ) * Y1(kl,ab)
rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,k,l) - ERI(p,q,l,k) &
+ 0d0*eh_Phi(p,q,k,l) - 0d0*eh_Phi(p,q,l,k) ) * Y1(kl,ab)
end do
end do
@ -131,8 +132,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
do c=nO+1,nOrb-nR
do d=c+1,nOrb-nR
cd = cd + 1
rho2(p,q,ij) = ( ERI(p,q,c,d) - ERI(p,q,d,c) ) * X2(cd,ij) &
+ ( eh_Phi(p,q,c,d) - eh_Phi(p,q,d,c) ) * X2(cd,ij)
rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,c,d) - ERI(p,q,d,c) &
+ 0d0*eh_Phi(p,q,c,d) - 0d0*eh_Phi(p,q,d,c) ) * X2(cd,ij)
end do
end do
@ -140,8 +141,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho2(p,q,ij) = ( ERI(p,q,k,l) - ERI(p,q,l,k) ) * Y2(kl,ij) &
+ ( eh_Phi(p,q,k,l) - eh_Phi(p,q,l,k) ) * Y2(kl,ij)
rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,k,l) - ERI(p,q,l,k) &
+ 0d0*eh_Phi(p,q,k,l) - 0d0*eh_Phi(p,q,l,k) ) * Y2(kl,ij)
end do
end do
end do

View File

@ -117,7 +117,7 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
allocate(old_eh_sing_Om(nS),old_eh_trip_Om(nS))
allocate(old_ee_sing_Om(nVVs),old_hh_sing_Om(nOOs))
allocate(old_ee_trip_Om(nVVt),old_hh_trip_Om(nOOt))
allocate(eh_sing_rho(nOrb,nOrb,nS+nS),eh_trip_rho(nOrb,nOrb,nS+nS))
allocate(eh_sing_rho(nOrb,nOrb,nS),eh_trip_rho(nOrb,nOrb,nS))
allocate(ee_sing_rho(nOrb,nOrb,nVVs),hh_sing_rho(nOrb,nOrb,nOOs))
allocate(ee_trip_rho(nOrb,nOrb,nVVt),hh_trip_rho(nOrb,nOrb,nOOt))
allocate(old_eh_sing_Phi(nOrb,nOrb,nOrb,nOrb),old_eh_trip_Phi(nOrb,nOrb,nOrb,nOrb))
@ -330,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)
@ -388,9 +388,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_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)
@ -439,7 +439,7 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
deallocate(eh_sing_rho,eh_trip_rho,ee_sing_rho,ee_trip_rho,hh_sing_rho,hh_trip_rho)
! TODO Once we will compute the blocks of kernel starting from the 4-tensors we can move the freeing up
! Memory allocation
allocate(eh_sing_rho(nOrb,nOrb,nS+nS),eh_trip_rho(nOrb,nOrb,nS+nS))
allocate(eh_sing_rho(nOrb,nOrb,nS),eh_trip_rho(nOrb,nOrb,nS))
allocate(ee_sing_rho(nOrb,nOrb,nVVs),hh_sing_rho(nOrb,nOrb,nOOs))
allocate(ee_trip_rho(nOrb,nOrb,nVVt),hh_trip_rho(nOrb,nOrb,nOOt))

View File

@ -7,7 +7,7 @@ subroutine R_eh_singlet_Phi(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_sing_Phi)
! 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)
double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS)
! Local variables
integer :: p,q,r,s
@ -24,11 +24,11 @@ subroutine R_eh_singlet_Phi(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_sing_Phi)
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
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)*eh_sing_rho(s,q,n)/eh_sing_Om(n)
end do
enddo
enddo

View File

@ -7,7 +7,7 @@ subroutine R_eh_triplet_Phi(nOrb,nC,nR,nS,eh_trip_Om,eh_trip_rho,eh_trip_Phi)
! 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)
double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS)
! Local variables
integer :: p,q,r,s
@ -24,11 +24,11 @@ subroutine R_eh_triplet_Phi(nOrb,nC,nR,nS,eh_trip_Om,eh_trip_rho,eh_trip_Phi)
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
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)*eh_trip_rho(s,q,n)/eh_trip_Om(n)
end do
enddo
enddo

View File

@ -35,7 +35,7 @@ subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nOOs,eh_sing_Phi,eh_trip_Phi,pp_sing_
kl = kl +1
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)
+ 0.5d0*eh_sing_Phi(i,j,l,k) - 1.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)))
@ -86,7 +86,7 @@ subroutine R_pp_singlet_Gamma_C(nOrb,nO,nR,nVVs,eh_sing_Phi,eh_trip_Phi,pp_sing_
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)
+ 0.5d0*eh_sing_Phi(a,b,d,c) - 1.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)))
@ -136,7 +136,7 @@ subroutine R_pp_singlet_Gamma_B(nOrb,nC,nO,nR,nOOs,nVVs,eh_sing_Phi,eh_trip_Phi,
ij = ij +1
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)
+ 0.5d0*eh_sing_Phi(a,b,j,i) - 1.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)))

View File

@ -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+nS)
double precision,intent(out) :: rho(nOrb,nOrb,nS)
rho(:,:,:) = 0d0
! !$OMP PARALLEL &
@ -27,24 +27,23 @@ subroutine R_eh_singlet_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
do ia=1,nS
jb = 0
do j=nC+1,nO
do b=nO+1,nOrb-nR
jb = jb + 1
X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb))
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
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
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
rho(p,q,ia) = rho(p,q,ia) + (2d0*ERI(q,j,p,b) - ERI(q,j,b,p) &
- 0d0*0.5d0*eh_sing_Phi(q,j,b,p) - 0d0*1.5d0*eh_trip_Phi(q,j,b,p) &
+ 0d0*0.5d0*pp_sing_Phi(q,j,p,b) + 0d0*1.5d0*pp_trip_Phi(q,j,p,b)) * X &
+ (2d0*ERI(q,b,p,j) - ERI(q,b,j,p) &
- 0d0*0.5d0*eh_sing_Phi(q,b,j,p) - 0d0*1.5d0*eh_trip_Phi(q,b,j,p) &
+ 0d0*0.5d0*pp_sing_Phi(q,b,p,j) + 0d0*1.5d0*pp_trip_Phi(q,b,p,j)) * Y
end do
@ -77,7 +76,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+nS)
double precision,intent(out) :: rho(nOrb,nOrb,nS)
rho(:,:,:) = 0d0
! !$OMP PARALLEL &
@ -87,28 +86,27 @@ 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
do ia=1,nS
jb = 0
do j=nC+1,nO
do b=nO+1,nOrb-nR
jb = jb + 1
X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb))
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
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
rho(p,q,ia) = rho(p,q,ia) + (- ERI(q,j,b,p) &
- 0d0*0.5d0*eh_sing_Phi(q,j,b,p) + 0d0*0.5d0*eh_trip_Phi(q,j,b,p) &
- 0d0*0.5d0*pp_sing_Phi(q,j,p,b) + 0d0*0.5d0*pp_trip_Phi(q,j,p,b)) * X &
+ (- ERI(q,b,j,p) &
- 0d0*0.5d0*eh_sing_Phi(q,b,j,p) + 0d0*0.5d0*eh_trip_Phi(q,b,j,p) &
- 0d0*0.5d0*pp_sing_Phi(q,b,p,j) + 0d0*0.5d0*pp_trip_Phi(q,b,p,j)) * Y
end do
end do
end do
end do
@ -175,9 +173,9 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
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))&
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) + ERI(p,q,d,c) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) - 0d0*1.5d0*eh_trip_Phi(p,q,c,d) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*1.5d0*eh_trip_Phi(p,q,d,c))&
*X1(cd,ab)/sqrt(1d0 + Kronecker_delta(c,d))
end do ! d
@ -189,9 +187,9 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
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))&
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) + ERI(p,q,l,k) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) - 0d0*1.5d0*eh_trip_Phi(p,q,k,l) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*1.5d0*eh_trip_Phi(p,q,l,k))&
*Y1(kl,ab)/sqrt(1d0 + Kronecker_delta(k,l))
end do ! l
end do ! k
@ -208,9 +206,9 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
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))&
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) + ERI(p,q,d,c) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) - 0d0*1.5d0*eh_trip_Phi(p,q,c,d) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*1.5d0*eh_trip_Phi(p,q,d,c))&
*X2(cd,ij)/sqrt(1d0 + Kronecker_delta(c,d))
end do ! d
end do ! c
@ -220,9 +218,9 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
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))&
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) + ERI(p,q,l,k) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) - 0d0*1.5d0*eh_trip_Phi(p,q,k,l) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*1.5d0*eh_trip_Phi(p,q,l,k))&
*Y2(kl,ij)/sqrt(1d0 + Kronecker_delta(k,l))
end do ! l
end do ! k
@ -292,9 +290,9 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
do d = c+1, 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) + 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)
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) - ERI(p,q,d,c) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) + 0d0*0.5d0*eh_trip_Phi(p,q,c,d) &
- 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*0.5d0*eh_trip_Phi(p,q,d,c) )*X1(cd,ab)
end do ! d
end do ! c
@ -305,9 +303,9 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
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) + 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)
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) - ERI(p,q,l,k) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) + 0d0*0.5d0*eh_trip_Phi(p,q,k,l) &
- 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*0.5d0*eh_trip_Phi(p,q,l,k) )*Y1(kl,ab)
end do ! l
end do ! k
end do ! b
@ -323,9 +321,9 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
do d = c+1, 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) + 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)
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) - ERI(p,q,d,c) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) + 0d0*0.5d0*eh_trip_Phi(p,q,c,d) &
- 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*0.5d0*eh_trip_Phi(p,q,d,c) )*X2(cd,ij)
end do ! d
end do ! c
@ -334,9 +332,9 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
do l = k+1, 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) + 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)
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) - ERI(p,q,l,k) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) + 0d0*0.5d0*eh_trip_Phi(p,q,k,l) &
- 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*0.5d0*eh_trip_Phi(p,q,l,k) )*Y2(kl,ij)
end do ! l
end do ! k