diff --git a/src/Parquet/README.md b/src/Parquet/README.md index 6a0fe5f..ee1ef1e 100644 --- a/src/Parquet/README.md +++ b/src/Parquet/README.md @@ -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$ diff --git a/src/Parquet/RParquet.f90 b/src/Parquet/RParquet.f90 index cdb1f07..b909a44 100644 --- a/src/Parquet/RParquet.f90 +++ b/src/Parquet/RParquet.f90 @@ -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 diff --git a/src/Parquet/R_eh_singlet_Gam.f90 b/src/Parquet/R_eh_singlet_Gam.f90 index 546bb5a..85e8cb0 100644 --- a/src/Parquet/R_eh_singlet_Gam.f90 +++ b/src/Parquet/R_eh_singlet_Gam.f90 @@ -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 diff --git a/src/Parquet/R_eh_triplet_Gam.f90 b/src/Parquet/R_eh_triplet_Gam.f90 index eb019c8..13b29d5 100644 --- a/src/Parquet/R_eh_triplet_Gam.f90 +++ b/src/Parquet/R_eh_triplet_Gam.f90 @@ -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 diff --git a/src/Parquet/R_pp_singlet_Gam.f90 b/src/Parquet/R_pp_singlet_Gam.f90 index 9bc0e7b..83e4826 100644 --- a/src/Parquet/R_pp_singlet_Gam.f90 +++ b/src/Parquet/R_pp_singlet_Gam.f90 @@ -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))) diff --git a/src/Parquet/R_pp_triplet_Gam.f90 b/src/Parquet/R_pp_triplet_Gam.f90 index 7c3efbb..9e12c3c 100644 --- a/src/Parquet/R_pp_triplet_Gam.f90 +++ b/src/Parquet/R_pp_triplet_Gam.f90 @@ -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 diff --git a/src/Parquet/R_screened_integrals.f90 b/src/Parquet/R_screened_integrals.f90 index 97c24c6..01b3f87 100644 --- a/src/Parquet/R_screened_integrals.f90 +++ b/src/Parquet/R_screened_integrals.f90 @@ -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