PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-04-01 06:21:37 +02:00
Code Releases Activity

update Parquet

Browse Source
This commit is contained in:
Antoine Marie 2025-03-20 16:27:04 +01:00
parent 4774ab082a
commit 4e7ce9dbb4
7 changed files with 318 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/Parquet/README.md
View File

@ -21,6 +21,7 @@ The hard-coded parameters are:
## TODO list
### Check
- [x] Initial ppRPA@HF eigenvalues checked with Ne DIP in Table 1 of ppBSE paper
- [ ] Comment m,s,t channels and perform ehBSE@$GW$ and ppBSE@$GW$
- [ ] Comment d,m channels and perform ehBSE@$GT$ and ppBSE@$GT$

99
src/Parquet/RParquet.f90
View File

@ -103,7 +103,7 @@ subroutine RParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body,nOrb,n
write(*,*)
write(*,*)'************ Solving initial linear-response problems ************'
write(*,*)'------------------------------------------------------------------'
!-------------------
! Density channel
!-------------------
@ -118,8 +118,7 @@ subroutine RParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body,nOrb,n
t = end_t - start_t
write(*,*)
write(*,'(A51,1X,F9.3,A8)') 'Total wall time for initial singlet phRPA problem =',t,' seconds'
!if (print_phLR) call print_excitation_energies('phRPA@RHF','singlet',nS,eh_sing_Om)
call print_excitation_energies('phRPA@RHF','singlet',nS,eh_sing_Om)
if (print_phLR) call print_excitation_energies('phRPA@RHF','singlet',nS,eh_sing_Om)
deallocate(Aph,Bph)
@ -137,8 +136,7 @@ subroutine RParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body,nOrb,n
t = end_t - start_t
write(*,*)
write(*,'(A51,1X,F9.3,A8)') 'Total wall time for initial triplet phRPA problem =',t,' seconds'
!if (print_phLR) call print_excitation_energies('phRPA@RHF','triplet',nS,eh_trip_Om)
call print_excitation_energies('phRPA@RHF','triplet',nS,eh_trip_Om)
if (print_phLR) call print_excitation_energies('phRPA@RHF','triplet',nS,eh_trip_Om)
deallocate(Aph,Bph)
@ -202,24 +200,50 @@ subroutine RParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body,nOrb,n
deallocate(eh_sing_Om,eh_trip_Om,ee_sing_Om,hh_sing_Om,ee_trip_Om,hh_trip_Om)
allocate(eh_sing_rho(nOrb,nOrb,nS))
call R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,sing_XpY,eh_sing_rho)
allocate(eh_sing_Gam(nOrb,nOrb,nOrb,nOrb))
eh_sing_Gam(:,:,:,:) = 0d0
call wall_time(start_t)
call R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,sing_XpY,eh_sing_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A51,1X,F9.3,A8)') 'Total wall time for building singlet eh screened integrals =',t,' seconds'
write(*,*)
deallocate(sing_XpY,sing_XmY)
deallocate(eh_sing_Gam)
allocate(eh_trip_rho(nOrb,nOrb,nS))
call R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,trip_XpY,eh_trip_rho)
allocate(eh_trip_Gam(nOrb,nOrb,nOrb,nOrb))
eh_trip_Gam(:,:,:,:) = 0d0
call wall_time(start_t)
call R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,trip_XpY,eh_trip_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A51,1X,F9.3,A8)') 'Total wall time for building triplet eh screened integrals =',t,' seconds'
write(*,*)
deallocate(trip_XpY,trip_XmY)
deallocate(eh_trip_Gam)
allocate(ee_sing_rho(nOrb,nOrb,nVVs),hh_sing_rho(nOrb,nOrb,nOOs))
allocate(pp_sing_Gam(nOrb,nOrb,nOrb,nOrb))
pp_sing_Gam(:,:,:,:) = 0d0
call wall_time(start_t)
call R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,pp_sing_Gam,X1s,Y1s,ee_sing_rho,X2s,Y2s,hh_sing_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A51,1X,F9.3,A8)') 'Total wall time for building singlet pp screened integrals =',t,' seconds'
write(*,*)
deallocate(X1s,Y1s,X2s,Y2s)
deallocate(pp_sing_Gam)
allocate(ee_trip_rho(nOrb,nOrb,nVVt),hh_trip_rho(nOrb,nOrb,nOOt))
allocate(pp_trip_Gam(nOrb,nOrb,nOrb,nOrb))
pp_trip_Gam(:,:,:,:) = 0d0
call wall_time(start_t)
call R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,pp_trip_Gam,X1t,Y1t,ee_trip_rho,X2t,Y2t,hh_trip_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A51,1X,F9.3,A8)') 'Total wall time for building triplet pp screened integrals =',t,' seconds'
write(*,*)
deallocate(X1t,Y1t,X2t,Y2t)
deallocate(pp_trip_Gam)
@ -398,6 +422,29 @@ subroutine RParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body,nOrb,n
old_ee_trip_Om(:) = ee_trip_Om(:)
old_hh_trip_Om(:) = hh_trip_Om(:)
!-------------------
! Compute effective interactions
!-------------------
allocate(eh_sing_Gam(nOrb,nOrb,nOrb,nOrb))
call wall_time(start_t)
call R_eh_singlet_Gamma(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, &
old_ee_sing_Om,ee_sing_rho,old_ee_trip_Om,ee_trip_rho, &
old_hh_sing_Om,hh_sing_rho,old_hh_trip_Om,hh_trip_rho, eh_sing_Gam)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A52,1X,F9.3,A8)') 'Total wall time for building eh singlet Gamma =',t,' seconds'
write(*,*)
allocate(eh_trip_Gam(nOrb,nOrb,nOrb,nOrb))
call wall_time(start_t)
call R_eh_triplet_Gamma(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, &
old_ee_sing_Om,ee_sing_rho,old_ee_trip_Om,ee_trip_rho, &
old_hh_sing_Om,hh_sing_rho,old_hh_trip_Om,hh_trip_rho, eh_trip_Gam)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A52,1X,F9.3,A8)') 'Total wall time for building eh triplet Gamma =',t,' seconds'
write(*,*)
allocate(pp_trip_Gam(nOrb,nOrb,nOrb,nOrb))
call wall_time(start_t)
call R_pp_triplet_Gamma(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho,pp_trip_Gam)
@ -412,24 +459,50 @@ subroutine RParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body,nOrb,n
t = end_t - start_t
write(*,'(A52,1X,F9.3,A8)') 'Total wall time for building pp singlet Gamma =',t,' seconds'
write(*,*)
!-------------------
! Deallocate
!-------------------
deallocate(eh_sing_Om,eh_trip_Om,ee_sing_Om,hh_sing_Om,ee_trip_Om,hh_trip_Om)
deallocate(eh_sing_rho,eh_trip_rho,ee_sing_rho,ee_trip_rho,hh_sing_rho,hh_trip_rho)
!-------------------
! Compute screened integrals
!-------------------
allocate(eh_sing_rho(nOrb,nOrb,nS))
call R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,sing_XpY,eh_sing_rho)
deallocate(sing_XpY,sing_XmY)
call wall_time(start_t)
call R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,sing_XpY,eh_sing_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A51,1X,F9.3,A8)') 'Total wall time for building singlet eh screened integrals =',t,' seconds'
write(*,*)
deallocate(sing_XpY,sing_XmY,eh_sing_Gam)
allocate(eh_trip_rho(nOrb,nOrb,nS))
call R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,trip_XpY,eh_trip_rho)
deallocate(trip_XpY,trip_XmY)
call wall_time(start_t)
call R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,trip_XpY,eh_trip_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A51,1X,F9.3,A8)') 'Total wall time for building triplet eh screened integrals =',t,' seconds'
write(*,*)
deallocate(trip_XpY,trip_XmY,eh_trip_Gam)
allocate(ee_sing_rho(nOrb,nOrb,nVVs),hh_sing_rho(nOrb,nOrb,nOOs))
call wall_time(start_t)
call R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,pp_sing_Gam,X1s,Y1s,ee_sing_rho,X2s,Y2s,hh_sing_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A51,1X,F9.3,A8)') 'Total wall time for building singlet pp screened integrals =',t,' seconds'
write(*,*)
deallocate(X1s,Y1s,X2s,Y2s,pp_sing_Gam)
allocate(ee_trip_rho(nOrb,nOrb,nVVt),hh_trip_rho(nOrb,nOrb,nOOt))
call wall_time(start_t)
call R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,pp_trip_Gam,X1t,Y1t,ee_trip_rho,X2t,Y2t,hh_trip_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A51,1X,F9.3,A8)') 'Total wall time for building triplet pp screened integrals =',t,' seconds'
write(*,*)
deallocate(X1t,Y1t,X2t,Y2t,pp_trip_Gam)
! Convergence criteria

81
src/Parquet/R_eh_singlet_Gam.f90
View File

@ -1,3 +1,76 @@
subroutine R_eh_singlet_Gamma(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, &
ee_sing_Om,ee_sing_rho,ee_trip_Om,ee_trip_rho, &
hh_sing_Om,hh_sing_rho,hh_trip_Om,hh_trip_rho, eh_sing_Gam)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt
double precision,intent(in) :: eh_sing_Om(nS)
double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS)
double precision,intent(in) :: eh_trip_Om(nS)
double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS)
double precision,intent(in) :: ee_sing_Om(nVVs)
double precision,intent(in) :: ee_sing_rho(nOrb,nOrb,nVVs)
double precision,intent(in) :: ee_trip_Om(nVVt)
double precision,intent(in) :: ee_trip_rho(nOrb,nOrb,nVVt)
double precision,intent(in) :: hh_sing_Om(nOOs)
double precision,intent(in) :: hh_sing_rho(nOrb,nOrb,nOOs)
double precision,intent(in) :: hh_trip_Om(nVVs)
double precision,intent(in) :: hh_trip_rho(nOrb,nOrb,nVVs)
! Local variables
integer :: p,q,r,s
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: eh_sing_Gam(nOrb,nOrb,nOrb,nOrb)
! Initialization
eh_sing_Gam(:,:,:,:) = 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_Gam(p,q,r,s) = eh_sing_Gam(p,q,r,s) &
+ eh_sing_rho(s,p,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(s,p,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n)
end do
do n=1,nVVs
eh_sing_Gam(p,q,r,s) = eh_sing_Gam(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,nOOs
eh_sing_Gam(p,q,r,s) = eh_sing_Gam(p,q,r,s) &
- hh_sing_rho(p,q,n)*hh_sing_rho(r,s,n)/hh_sing_Om(n)
end do
do n=1,nVVt
eh_sing_Gam(p,q,r,s) = eh_sing_Gam(p,q,r,s) &
+ 3d0 * ee_trip_rho(p,q,n)*ee_trip_rho(r,s,n)/ee_trip_Om(n)
end do
do n=1,nOOt
eh_sing_Gam(p,q,r,s) = eh_sing_Gam(p,q,r,s) &
- 3d0 * 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_eh_singlet_Gamma
subroutine R_eh_singlet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, &
ee_sing_Om,ee_sing_rho,ee_trip_Om,ee_trip_rho, &
@ -45,8 +118,8 @@ subroutine R_eh_singlet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
do n=1,nS
eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) &
+ eh_sing_rho(a,b,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(a,b,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n)
+ eh_sing_rho(b,a,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(b,a,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n)
end do
do n=1,nVVs
@ -123,8 +196,8 @@ subroutine R_eh_singlet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
do n=1,nS
eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) &
+ eh_sing_rho(a,j,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(a,j,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n)
+ eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n)
end do
do n=1,nVVs

81
src/Parquet/R_eh_triplet_Gam.f90
View File

@ -1,3 +1,76 @@
subroutine R_eh_triplet_Gamma(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, &
ee_sing_Om,ee_sing_rho,ee_trip_Om,ee_trip_rho, &
hh_sing_Om,hh_sing_rho,hh_trip_Om,hh_trip_rho, eh_trip_Gam)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt
double precision,intent(in) :: eh_sing_Om(nS)
double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS)
double precision,intent(in) :: eh_trip_Om(nS)
double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS)
double precision,intent(in) :: ee_sing_Om(nVVs)
double precision,intent(in) :: ee_sing_rho(nOrb,nOrb,nVVs)
double precision,intent(in) :: ee_trip_Om(nVVt)
double precision,intent(in) :: ee_trip_rho(nOrb,nOrb,nVVt)
double precision,intent(in) :: hh_sing_Om(nOOs)
double precision,intent(in) :: hh_sing_rho(nOrb,nOrb,nOOs)
double precision,intent(in) :: hh_trip_Om(nVVs)
double precision,intent(in) :: hh_trip_rho(nOrb,nOrb,nVVs)
! Local variables
integer :: p,q,r,s
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: eh_trip_Gam(nOrb,nOrb,nOrb,nOrb)
! Initialization
eh_trip_Gam(:,:,:,:) = 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_Gam(p,q,r,s) = eh_trip_Gam(p,q,r,s) &
+ eh_sing_rho(s,p,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) &
- eh_trip_rho(s,p,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n)
end do
do n=1,nVVs
eh_trip_Gam(p,q,r,s) = eh_trip_Gam(p,q,r,s) &
- 0d0*ee_sing_rho(p,q,n) * ee_sing_rho(r,s,n)/ee_sing_Om(n)
end do
do n=1,nOOs
eh_trip_Gam(p,q,r,s) = eh_trip_Gam(p,q,r,s) &
+ 0d0*hh_sing_rho(p,q,n) * hh_sing_rho(r,s,n)/hh_sing_Om(n)
end do
do n=1,nVVt
eh_trip_Gam(p,q,r,s) = eh_trip_Gam(p,q,r,s) &
+ 0d0*ee_trip_rho(p,q,n) * ee_trip_rho(r,s,n)/ee_trip_Om(n)
end do
do n=1,nOOt
eh_trip_Gam(p,q,r,s) = eh_trip_Gam(p,q,r,s) &
- 0d0*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_eh_triplet_Gamma
subroutine R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, &
ee_sing_Om,ee_sing_rho,ee_trip_Om,ee_trip_rho, &
@ -45,8 +118,8 @@ subroutine R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
do n=1,nS
eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) &
+ eh_sing_rho(a,b,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) &
- eh_trip_rho(a,b,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n)
+ eh_sing_rho(b,a,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) &
- eh_trip_rho(b,a,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n)
end do
do n=1,nVVs
@ -123,8 +196,8 @@ subroutine R_eh_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
do n=1,nS
eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) &
+ eh_sing_rho(a,j,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
- eh_trip_rho(a,j,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n)
+ eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
- eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n)
end do
do n=1,nVVs

32
src/Parquet/R_pp_singlet_Gam.f90
View File

@ -33,10 +33,10 @@ subroutine R_pp_singlet_Gamma(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh
do n=1,nS
pp_sing_Gam(p,q,r,s) = pp_sing_Gam(p,q,r,s) &
- eh_sing_rho(p,r,n)*eh_sing_rho(q,s,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(p,r,n)*eh_trip_rho(q,s,n)/eh_trip_Om(n) &
- eh_sing_rho(p,s,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(p,s,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n)
- eh_sing_rho(r,p,n)*eh_sing_rho(q,s,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(r,p,n)*eh_trip_rho(q,s,n)/eh_trip_Om(n) &
- eh_sing_rho(s,p,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(s,p,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n)
end do
pp_sing_Gam(p,q,r,s) = pp_sing_Gam(p,q,r,s)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(r,s)))
@ -91,10 +91,10 @@ subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOs,eh_sing_Om,eh_sing_rho,
do n=1,nS
pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl) &
- eh_sing_rho(i,k,n)*eh_sing_rho(j,l,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(i,k,n)*eh_trip_rho(j,l,n)/eh_trip_Om(n) &
- eh_sing_rho(i,l,n)*eh_sing_rho(j,k,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(i,l,n)*eh_trip_rho(j,k,n)/eh_trip_Om(n)
- eh_sing_rho(k,i,n)*eh_sing_rho(j,l,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(k,i,n)*eh_trip_rho(j,l,n)/eh_trip_Om(n) &
- eh_sing_rho(l,i,n)*eh_sing_rho(j,k,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(l,i,n)*eh_trip_rho(j,k,n)/eh_trip_Om(n)
end do
pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
@ -149,10 +149,10 @@ subroutine R_pp_singlet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVs,eh_sing_Om,eh_sing_rho,
do n=1,nS
pp_sing_Gam_C(ab,cd) = pp_sing_Gam_C(ab,cd) &
- eh_sing_rho(a,c,n)*eh_sing_rho(b,d,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(a,c,n)*eh_trip_rho(b,d,n)/eh_trip_Om(n) &
- eh_sing_rho(a,d,n)*eh_sing_rho(b,c,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(a,d,n)*eh_trip_rho(b,c,n)/eh_trip_Om(n)
- eh_sing_rho(c,a,n)*eh_sing_rho(b,d,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(c,a,n)*eh_trip_rho(b,d,n)/eh_trip_Om(n) &
- eh_sing_rho(c,a,n)*eh_sing_rho(b,c,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(d,a,n)*eh_trip_rho(b,c,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)))
@ -207,10 +207,10 @@ subroutine R_pp_singlet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,eh_sing_Om,eh_sing
do n=1,nS
pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij) &
- eh_sing_rho(a,i,n)*eh_sing_rho(b,j,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(a,i,n)*eh_trip_rho(b,j,n)/eh_trip_Om(n) &
- eh_sing_rho(a,j,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(a,j,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n)
- eh_sing_rho(i,a,n)*eh_sing_rho(b,j,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(i,a,n)*eh_trip_rho(b,j,n)/eh_trip_Om(n) &
- eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
+ 3d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n)
end do
pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))

32
src/Parquet/R_pp_triplet_Gam.f90
View File

@ -33,10 +33,10 @@ subroutine R_pp_triplet_Gamma(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh
do n=1,nS
pp_trip_Gam(p,q,r,s) = pp_trip_Gam(p,q,r,s) &
- eh_sing_rho(p,r,n)*eh_sing_rho(q,s,n)/eh_sing_Om(n) &
- eh_trip_rho(p,r,n)*eh_trip_rho(q,s,n)/eh_trip_Om(n) &
+ eh_sing_rho(p,s,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) &
+ eh_trip_rho(p,s,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n)
- eh_sing_rho(r,p,n)*eh_sing_rho(q,s,n)/eh_sing_Om(n) &
- eh_trip_rho(r,p,n)*eh_trip_rho(q,s,n)/eh_trip_Om(n) &
+ eh_sing_rho(s,p,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) &
+ eh_trip_rho(s,p,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n)
end do
end do
@ -89,10 +89,10 @@ subroutine R_pp_triplet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOt,eh_sing_Om,eh_sing_rho,
do n=1,nS
pp_trip_Gam_D(ij,kl) = pp_trip_Gam_D(ij,kl) &
- eh_sing_rho(i,k,n)*eh_sing_rho(j,l,n)/eh_sing_Om(n) &
- eh_trip_rho(i,k,n)*eh_trip_rho(j,l,n)/eh_trip_Om(n) &
+ eh_sing_rho(i,l,n)*eh_sing_rho(j,k,n)/eh_sing_Om(n) &
+ eh_trip_rho(i,l,n)*eh_trip_rho(j,k,n)/eh_trip_Om(n)
- eh_sing_rho(k,i,n)*eh_sing_rho(j,l,n)/eh_sing_Om(n) &
- eh_trip_rho(k,i,n)*eh_trip_rho(j,l,n)/eh_trip_Om(n) &
+ eh_sing_rho(l,i,n)*eh_sing_rho(j,k,n)/eh_sing_Om(n) &
+ eh_trip_rho(l,i,n)*eh_trip_rho(j,k,n)/eh_trip_Om(n)
end do
end do
@ -145,10 +145,10 @@ subroutine R_pp_triplet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVt,eh_sing_Om,eh_sing_rho,
do n=1,nS
pp_trip_Gam_C(ab,cd) = pp_trip_Gam_C(ab,cd) &
- eh_sing_rho(a,c,n)*eh_sing_rho(b,d,n)/eh_sing_Om(n) &
- eh_trip_rho(a,c,n)*eh_trip_rho(b,d,n)/eh_trip_Om(n) &
+ eh_sing_rho(a,d,n)*eh_sing_rho(b,c,n)/eh_sing_Om(n) &
+ eh_trip_rho(a,d,n)*eh_trip_rho(b,c,n)/eh_trip_Om(n)
- eh_sing_rho(c,a,n)*eh_sing_rho(b,d,n)/eh_sing_Om(n) &
- eh_trip_rho(c,a,n)*eh_trip_rho(b,d,n)/eh_trip_Om(n) &
+ eh_sing_rho(d,a,n)*eh_sing_rho(b,c,n)/eh_sing_Om(n) &
+ eh_trip_rho(d,a,n)*eh_trip_rho(b,c,n)/eh_trip_Om(n)
end do
end do
@ -201,10 +201,10 @@ subroutine R_pp_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOt,nVVt,eh_sing_Om,eh_sing
do n=1,nS
pp_trip_Gam_B(ab,ij) = pp_trip_Gam_B(ab,ij) &
- eh_sing_rho(a,i,n)*eh_sing_rho(b,j,n)/eh_sing_Om(n) &
- eh_trip_rho(a,i,n)*eh_trip_rho(b,j,n)/eh_trip_Om(n) &
+ eh_sing_rho(a,j,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
+ eh_trip_rho(a,j,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n)
- eh_sing_rho(i,a,n)*eh_sing_rho(b,j,n)/eh_sing_Om(n) &
- eh_trip_rho(i,a,n)*eh_trip_rho(b,j,n)/eh_trip_Om(n) &
+ eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
+ eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n)
end do
end do

88
src/Parquet/R_screened_integrals.f90
View File

@ -1,4 +1,4 @@
subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,XpY,rho)
! Compute excitation densities
implicit none
@ -6,6 +6,7 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
! Input variables
integer,intent(in) :: nOrb,nC,nO,nR,nS
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_sing_Gam(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: XpY(nS,nS)
! Local variables
@ -15,11 +16,11 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
double precision,intent(out) :: rho(nOrb,nOrb,nS)
rho(:,:,:) = 0d0
!$OMP PARALLEL &
!$OMP SHARED(nC,nOrb,nR,nO,nS,rho,ERI,XpY) &
!$OMP PRIVATE(q,p,jb,ia) &
!$OMP DEFAULT(NONE)
!$OMP DO
! !$OMP PARALLEL &
! !$OMP SHARED(nC,nOrb,nR,nO,nS,rho,ERI,XpY,eh_sing_Gam) &
! !$OMP PRIVATE(q,p,jb,ia) &
! !$OMP DEFAULT(NONE)
! !$OMP DO
do q=nC+1,nOrb-nR
do p=nC+1,nOrb-nR
jb = 0
@ -27,18 +28,19 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
do b=nO+1,nOrb-nR
jb = jb + 1
do ia=1,nS
rho(p,q,ia) = rho(p,q,ia) + (2d0*ERI(p,j,q,b) - ERI(p,j,b,q))*XpY(ia,jb)
rho(p,q,ia) = rho(p,q,ia) + (2d0*ERI(p,j,q,b) - ERI(p,j,b,q))*XpY(ia,jb) &
+ 0d0*eh_sing_Gam(p,j,q,b) * XpY(ia,jb)
end do
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine R_eh_singlet_screened_integral
subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,XpY,rho)
! Compute excitation densities
implicit none
@ -46,6 +48,7 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
! Input variables
integer,intent(in) :: nOrb,nC,nO,nR,nS
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Gam(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: XpY(nS,nS)
! Local variables
@ -55,11 +58,11 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
double precision,intent(out) :: rho(nOrb,nOrb,nS)
rho(:,:,:) = 0d0
!$OMP PARALLEL &
!$OMP SHARED(nC,nOrb,nR,nO,nS,rho,ERI,XpY) &
!$OMP PRIVATE(q,p,jb,ia) &
!$OMP DEFAULT(NONE)
!$OMP DO
! !$OMP PARALLEL &
! !$OMP SHARED(nC,nOrb,nR,nO,nS,rho,ERI,XpY,eh_trip_Gam) &
! !$OMP PRIVATE(q,p,jb,ia) &
! !$OMP DEFAULT(NONE)
! !$OMP DO
do q=nC+1,nOrb-nR
do p=nC+1,nOrb-nR
jb = 0
@ -67,14 +70,14 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
do b=nO+1,nOrb-nR
jb = jb + 1
do ia=1,nS
rho(p,q,ia) = rho(p,q,ia) - ERI(p,j,b,q)*XpY(ia,jb)
rho(p,q,ia) = rho(p,q,ia) - ERI(p,j,b,q)*XpY(ia,jb) + 0d0*eh_trip_Gam(p,j,q,b) * XpY(ia,jb)
end do
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine R_eh_triplet_screened_integral
@ -118,10 +121,10 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_sing_G
rho1(:,:,:) = 0d0
rho2(:,:,:) = 0d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
!$OMP SHARED(nC, nOrb, nR, nO, rho1, rho2, ERI, pp_sing_Gam, X1, Y1, X2, Y2)
!$OMP DO COLLAPSE(2)
! !$OMP PARALLEL DEFAULT(NONE) &
! !$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
! !$OMP SHARED(nC, nOrb, nR, nO, rho1, rho2, ERI, pp_sing_Gam, X1, Y1, X2, Y2)
! !$OMP DO COLLAPSE(2)
do q=nC+1,nOrb-nR
do p=nC+1,nOrb-nR
@ -134,10 +137,9 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_sing_G
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) + 0d0*pp_sing_Gam(p,q,c,d))*X1(cd,ab)/ &
sqrt(1d0 + Kronecker_delta(c,d))
+ (ERI(p,q,c,d) + ERI(p,q,d,c) + pp_sing_Gam(p,q,c,d))*X1(cd,ab) &
/sqrt(1d0 + Kronecker_delta(c,d))
end do
end do
@ -146,8 +148,8 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_sing_G
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) + 0d0*pp_sing_Gam(p,q,k,l))*Y1(kl,ab)/ &
sqrt(1d0 + Kronecker_delta(k,l))
+ (ERI(p,q,k,l) + ERI(p,q,l,k) + pp_sing_Gam(p,q,k,l))*Y1(kl,ab) &
/sqrt(1d0 + Kronecker_delta(k,l))
end do
end do
@ -163,8 +165,8 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_sing_G
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) + 0d0*pp_sing_Gam(p,q,c,d))*X2(cd,ij)/ &
sqrt(1d0 + Kronecker_delta(c,d))
+ (ERI(p,q,c,d) + ERI(p,q,d,c) + pp_sing_Gam(p,q,c,d))*X2(cd,ij) &
/sqrt(1d0 + Kronecker_delta(c,d))
end do
end do
@ -173,16 +175,16 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_sing_G
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) + 0d0*pp_sing_Gam(p,q,k,l))*Y2(kl,ij)/ &
sqrt(1d0 + Kronecker_delta(k,l))
+ (ERI(p,q,k,l) + ERI(p,q,l,k) + 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
end do
!$OMP END DO
!$OMP END PARALLEL
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine R_pp_singlet_screened_integral
@ -224,10 +226,10 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_trip_G
dim_1 = (nOrb - nO) * (nOrb - nO - 1) / 2
dim_2 = nO * (nO - 1) / 2
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
!$OMP SHARED(nC, nOrb, nR, nO, rho1, rho2, ERI, pp_trip_Gam, X1, Y1, X2, Y2)
!$OMP DO COLLAPSE(2)
! !$OMP PARALLEL DEFAULT(NONE) &
! !$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
! !$OMP SHARED(nC, nOrb, nR, nO, rho1, rho2, ERI, pp_trip_Gam, X1, Y1, X2, Y2)
! !$OMP DO COLLAPSE(2)
do q = nC+1, nOrb-nR
do p = nC+1, nOrb-nR
ab = 0
@ -242,7 +244,7 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_trip_G
cd = cd + 1
rho1(p,q,ab) = rho1(p,q,ab) &
+ (ERI(p,q,c,d) - ERI(p,q,d,c) + 0d0*pp_trip_Gam(p,q,c,d))*X1(cd,ab)
+ (ERI(p,q,c,d) - ERI(p,q,d,c) + pp_trip_Gam(p,q,c,d))*X1(cd,ab)
end do ! d
end do ! c
@ -253,7 +255,7 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_trip_G
kl = kl + 1
rho1(p,q,ab) = rho1(p,q,ab) &
+ (ERI(p,q,k,l) - ERI(p,q,l,k) + 0d0*pp_trip_Gam(p,q,k,l))*Y1(kl,ab)
+ (ERI(p,q,k,l) - ERI(p,q,l,k) + pp_trip_Gam(p,q,k,l))*Y1(kl,ab)
end do ! l
end do ! k
end do ! b
@ -270,7 +272,7 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_trip_G
cd = cd + 1
rho2(p,q,ij) = rho2(p,q,ij) &
+ (ERI(p,q,c,d) - ERI(p,q,d,c) + 0d0*pp_trip_Gam(p,q,c,d))*X2(cd,ij)
+ (ERI(p,q,c,d) - ERI(p,q,d,c) + pp_trip_Gam(p,q,c,d))*X2(cd,ij)
end do ! d
end do ! c
@ -280,14 +282,14 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_trip_G
kl = kl + 1
rho2(p,q,ij) = rho2(p,q,ij) &
+ (ERI(p,q,k,l) - ERI(p,q,l,k) + 0d0*pp_trip_Gam(p,q,k,l))*Y2(kl,ij)
+ (ERI(p,q,k,l) - ERI(p,q,l,k) + pp_trip_Gam(p,q,k,l))*Y2(kl,ij)
end do ! l
end do ! k
end do ! j
end do ! i
end do ! p
end do ! q
!$OMP END DO
!$OMP END PARALLEL
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine R_pp_triplet_screened_integral