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

spin adaptation checked up to screened integrals, ie iteration 3

Browse Source
This commit is contained in:
Antoine Marie 2025-03-24 21:21:44 +01:00
parent 36685fcc73
commit e0fe498b51
16 changed files with 866 additions and 94 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/GT/RGTeh_SigC.f90
View File

@ -35,7 +35,7 @@ double precision function RGTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
do i=nC+1,nO
do m=1,nS
eps = w - e(i) + Om(m)
num = rhoL(i,p,m)*rhoR(i,p,m)
num = rhoL(i,p,m)*rhoL(i,p,m)
RGTeh_SigC = RGTeh_SigC + num*eps/(eps**2 + eta**2)
end do
end do
@ -45,7 +45,7 @@ double precision function RGTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
do a=nO+1,nBas-nR
do m=1,nS
eps = w - e(a) - Om(m)
num = rhoL(p,a,m)*rhoR(p,a,m)
num = rhoL(p,a,m)*rhoL(p,a,m)
RGTeh_SigC = RGTeh_SigC + num*eps/(eps**2 + eta**2)
end do
end do

4
src/GT/RGTeh_dSigC.f90
View File

@ -35,7 +35,7 @@ double precision function RGTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR
do i=nC+1,nO
do m=1,nS
eps = w - e(i) + Om(m)
num = rhoL(i,p,m)*rhoR(i,p,m)
num = rhoL(i,p,m)*rhoL(i,p,m)
RGTeh_dSigC = RGTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
@ -45,7 +45,7 @@ double precision function RGTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR
do a=nO+1,nBas-nR
do m=1,nS
eps = w - e(a) - Om(m)
num = rhoL(p,a,m)*rhoR(p,a,m)
num = rhoL(p,a,m)*rhoL(p,a,m)
RGTeh_dSigC = RGTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do

4
src/GT/RGTeh_self_energy_diag.f90
View File

@ -46,7 +46,7 @@ subroutine RGTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Si
do m=1,nS
eps = e(p) - e(i) + Om(m)
num = rhoL(i,p,m)*rhoR(i,p,m)
num = rhoL(i,p,m)*rhoL(i,p,m)
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
@ -61,7 +61,7 @@ subroutine RGTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Si
do m=1,nS
eps = e(p) - e(a) - Om(m)
num = rhoL(p,a,m)*rhoR(p,a,m)
num = rhoL(p,a,m)*rhoL(p,a,m)
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2

2
src/LR/phGLR_A.f90
View File

@ -52,5 +52,5 @@ subroutine phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
end do
end do
end do
end subroutine

178
src/Parquet/GParquet.f90
View File

@ -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
@ -22,8 +22,6 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
! Local variables
integer :: ispin
integer :: n_it_1b,n_it_2b
double precision :: err_1b,err_2b
double precision :: err_eh, err_hh, err_ee
@ -155,9 +153,9 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
write(*,*) 'Computing eh effective interaction...'
call wall_time(start_t)
!call R_eh_Gamma(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
!call R_eh_Gamma(nOrb,nC,nO,nV,nR,nS,nOO,nVV, &
! old_eh_Om,eh_rho,old_ee_Om,ee_rho,old_hh_Om,hh_rho, &
! eh_trip_Gam)
! eh_Gam)
call wall_time(end_t)
t = end_t - start_t
@ -168,14 +166,180 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
write(*,*) 'Computing pp effective interaction...'
call wall_time(start_t)
!call R_pp_Gamma(nOrb,nC,nR,nS,old_eh_Om,eh_rho,pp_Gam)
call G_pp_Gamma(nOrb,nC,nO,nV,nR,nS,old_eh_Om,eh_rho,pp_Gam)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A50,1X,F9.3,A8)') 'Wall time for pp Gamma =',t,' seconds'
write(*,*)
!-----------------!
! eh channel !
!-----------------!
write(*,*)' ------------------------------'
write(*,*)' | Diagonalizing ehBSE |'
write(*,*)' ------------------------------'
write(*,*)
allocate(Aph(nS,nS),Bph(nS,nS),eh_Om(nS),XpY(nS,nS),XmY(nS,nS),eh_Gam_A(nS,nS),eh_Gam_B(nS,nS))
Aph(:,:) = 0d0
Bph(:,:) = 0d0
call wall_time(start_t)
call phGLR_A(.false.,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
if(.not.TDA) call phGLR_B(.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
if(n_it_2b == 1) then
eh_Gam_A(:,:) = 0d0
eh_Gam_B(:,:) = 0d0
else
call G_eh_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOO,nVV, &
old_eh_Om,eh_rho,old_ee_Om,ee_rho,old_hh_Om,hh_rho, &
eh_Gam_A)
if(.not.TDA) call G_eh_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOO,nVV, &
old_eh_Om,eh_rho,old_ee_Om,ee_rho,old_hh_Om,hh_rho, &
eh_Gam_B)
end if
Aph(:,:) = Aph(:,:) + eh_Gam_A(:,:)
Bph(:,:) = Bph(:,:) + eh_Gam_B(:,:)
call phGLR(TDA,nS,Aph,Bph,EcRPA,eh_Om,XpY,XmY)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A50,1X,F9.3,A8)') 'Wall time for phBSE =',t,' seconds'
write(*,*)
if(print_phLR) call print_excitation_energies('phBSE@Parquet','eh generalized',nS,eh_Om)
err_eh = maxval(abs(old_eh_Om - eh_Om))
deallocate(Aph,Bph,eh_Gam_A,eh_Gam_B)
!-----------------!
! pp channel !
!-----------------!
write(*,*)' ------------------------------'
write(*,*)' | Diagonalizing ppBSE |'
write(*,*)' ------------------------------'
write(*,*)
allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), &
ee_Om(nVV),X1(nVV,nVV),Y1(nOO,nVV), &
hh_Om(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
pp_Gam_B(nVV,nOO),pp_Gam_C(nVV,nVV),pp_Gam_D(nOO,nOO))
Bpp(:,:) = 0d0
Cpp(:,:) = 0d0
Dpp(:,:) = 0d0
call wall_time(start_t)
if(.not.TDA) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp)
call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp)
if(n_it_2b == 1) then
pp_Gam_B(:,:) = 0d0
pp_Gam_C(:,:) = 0d0
pp_Gam_D(:,:) = 0d0
else
if(.not.TDA) call G_pp_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOO,nVV,old_eh_Om,eh_rho,pp_Gam_B)
call G_pp_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVV,old_eh_Om,eh_rho,pp_Gam_C)
call G_pp_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOO,old_eh_Om,eh_rho,pp_Gam_D)
end if
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)
t = end_t - start_t
write(*,'(A50,1X,F9.3,A8)') 'Wall time for ppBSE =',t,' seconds'
write(*,*)
if(print_ppLR) call print_excitation_energies('ppBSE@Parquet','2p generalized',nVV,ee_Om)
if(print_ppLR) call print_excitation_energies('ppBSE@Parquet','2h generalized',nOO,hh_Om)
err_ee = maxval(abs(old_ee_Om - ee_Om))
err_hh = maxval(abs(old_hh_Om - hh_Om))
deallocate(Bpp,Cpp,Dpp,pp_Gam_B,pp_Gam_C,pp_Gam_D)
write(*,*) '----------------------------------------'
write(*,*) ' Two-body convergence '
write(*,*) '----------------------------------------'
write(*,'(1X,A30,F10.6)')'Error for eh channel = ',err_eh
write(*,'(1X,A30,F10.6)')'Error for pp channel = ',max(err_ee,err_hh)
write(*,*) '----------------------------------------'
write(*,*)
!----------!
! Updating !
!----------!
old_eh_Om(:) = eh_Om(:)
old_ee_Om(:) = ee_Om(:)
old_hh_Om(:) = hh_Om(:)
deallocate(eh_Om,ee_Om,hh_Om)
!----------------------------!
! Compute screened integrals !
!----------------------------!
! Free memory
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))
allocate(ee_rho(nOrb,nOrb,nVV),hh_rho(nOrb,nOrb,nOO))
! Build singlet eh integrals
write(*,*) 'Computing eh screened integrals...'
call wall_time(start_t)
call G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Gam,XpY,XmY,eh_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for eh integrals =',t,' seconds'
write(*,*)
! Done with eigenvectors and kernel
deallocate(XpY,XmY,eh_Gam)
! Build singlet pp integrals
write(*,*) 'Computing pp screened integrals...'
call wall_time(start_t)
call G_pp_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_Gam,X1,Y1,ee_rho,X2,Y2,hh_rho)
call wall_time(end_t)
t = end_t - start_t
! Done with eigenvectors and kernel
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for pp integrals =',t,' seconds'
write(*,*)
deallocate(X1,Y1,X2,Y2,pp_Gam)
! Convergence criteria
err_2b = max(err_eh,err_ee,err_hh)
call wall_time(end_2b)
t_2b = end_2b - start_2b

186
src/Parquet/G_eh_Gam.f90 Normal file
View File

@ -0,0 +1,186 @@
subroutine G_eh_Gamma(nOrb,nC,nO,nV,nR,nS,nOO,nVV, &
eh_Om,eh_rho,ee_Om,ee_rho,hh_Om,hh_rho, &
eh_Gam)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOO,nVV
double precision,intent(in) :: eh_Om(nS)
double precision,intent(in) :: eh_rho(nOrb,nOrb,nS)
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
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: eh_Gam(nOrb,nOrb,nOrb,nOrb)
! Initialization
eh_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 G_eh_Gamma
subroutine G_eh_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOO,nVV, &
eh_Om,eh_rho,ee_Om,ee_rho,hh_Om,hh_rho, &
eh_Gam_A)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOO,nVV
double precision,intent(in) :: eh_Om(nS)
double precision,intent(in) :: eh_rho(nOrb,nOrb,nS)
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 :: i,a,j,b
integer :: ia,jb
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: eh_Gam_A(nS,nS)
! Initialization
eh_Gam_A(:,:) = 0d0
ia = 0
do i=nC+1,nO
do a=nO+1,nOrb-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,norb-nR
jb = jb + 1
do n=1,nS
eh_Gam_A(ia,jb) = eh_Gam_A(ia,jb) &
+ eh_rho(b,a,n)*eh_rho(j,i,n)/eh_Om(n) &
+ eh_rho(a,b,n)*eh_rho(i,j,n)/eh_Om(n)
end do
do n=1,nVV
eh_Gam_A(ia,jb) = eh_Gam_A(ia,jb) &
+ 2d0 * ee_rho(a,j,n)*ee_rho(i,b,n)/ee_Om(n)
end do
do n=1,nOO
eh_Gam_A(ia,jb) = eh_Gam_A(ia,jb) &
- 2d0 * hh_rho(a,j,n)*hh_rho(i,b,n)/hh_Om(n)
end do
enddo
enddo
enddo
enddo
end subroutine G_eh_Gamma_A
subroutine G_eh_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOO,nVV, &
eh_Om,eh_rho,ee_Om,ee_rho,hh_Om,hh_rho, &
eh_Gam_B)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOO,nVV
double precision,intent(in) :: eh_Om(nS)
double precision,intent(in) :: eh_rho(nOrb,nOrb,nS)
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 :: i,a,j,b
integer :: ia,jb
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: eh_Gam_B(nS,nS)
! Initialization
eh_Gam_B(:,:) = 0d0
ia = 0
do i=nC+1,nO
do a=nO+1,nOrb-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,norb-nR
jb = jb + 1
do n=1,nS
eh_Gam_B(ia,jb) = eh_Gam_B(ia,jb) &
+ eh_rho(j,a,n)*eh_rho(b,i,n)/eh_Om(n) &
+ eh_rho(a,j,n)*eh_rho(i,b,n)/eh_Om(n)
end do
do n=1,nVV
eh_Gam_B(ia,jb) = eh_Gam_B(ia,jb) &
+ 2d0 * ee_rho(a,b,n)*ee_rho(i,j,n)/ee_Om(n)
end do
do n=1,nOO
eh_Gam_B(ia,jb) = eh_Gam_B(ia,jb) &
- 2d0 * hh_rho(a,b,n)*hh_rho(i,j,n)/hh_Om(n)
end do
enddo
enddo
enddo
enddo
end subroutine G_eh_Gamma_B

211
src/Parquet/G_pp_Gam.f90 Normal file
View File

@ -0,0 +1,211 @@
subroutine G_pp_Gamma(nOrb,nC,nO,nV,nR,nS,eh_Om,eh_rho,pp_Gam)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS
double precision,intent(in) :: eh_Om(nS)
double precision,intent(in) :: eh_rho(nOrb,nOrb,nS)
! Local variables
integer :: p,q,r,s
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_Gam(nOrb,nOrb,nOrb,nOrb)
! Initialization
pp_Gam(:,:,:,:) = 0d0
! !$OMP PARALLEL DEFAULT(NONE) &
! !$OMP PRIVATE(a, b, ab, i, j, ij, n) &
! !$OMP SHARED(nC, nOrb, nO, nS, pp_trip_Gam_B, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om)
! !$OMP DO COLLAPSE(2)
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
pp_Gam(p,q,r,s) = pp_Gam(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) &
+ eh_rho(s,p,n)*eh_rho(q,r,n)/eh_Om(n) &
+ eh_rho(p,s,n)*eh_rho(r,q,n)/eh_Om(n)
end do
end do
end do
end do
end do
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine
subroutine G_pp_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOO,eh_Om,eh_rho,pp_Gam_D)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOO
double precision,intent(in) :: eh_Om(nS)
double precision,intent(in) :: eh_rho(nOrb,nOrb,nS)
! Local variables
integer :: i,j,k,l
integer :: ij,kl
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_Gam_D(nOO,nOO)
! Initialization
pp_Gam_D(:,:) = 0d0
! !$OMP PARALLEL DEFAULT(NONE) &
! !$OMP PRIVATE(a, b, ab, i, j, ij, n) &
! !$OMP SHARED(nC, nOrb, nO, nS, pp_trip_Gam_B, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om)
! !$OMP DO COLLAPSE(2)
ij = 0
do i=nC+1,nO
do j=i+1,nO
ij = ij + 1
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl +1
do n=1,nS
pp_Gam_D(ij,kl) = pp_Gam_D(ij,kl) &
- eh_rho(k,i,n)*eh_rho(j,l,n)/eh_Om(n) &
- eh_rho(i,k,n)*eh_rho(l,j,n)/eh_Om(n) &
+ eh_rho(l,i,n)*eh_rho(j,k,n)/eh_Om(n) &
+ eh_rho(i,l,n)*eh_rho(k,j,n)/eh_Om(n)
end do
end do
end do
end do
end do
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine
subroutine G_pp_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVV,eh_Om,eh_rho,pp_Gam_C)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nVV
double precision,intent(in) :: eh_Om(nS)
double precision,intent(in) :: eh_rho(nOrb,nOrb,nS)
! Local variables
integer :: a,b,c,d
integer :: ab,cd
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_Gam_C(nVV,nVV)
! Initialization
pp_Gam_C(:,:) = 0d0
! !$OMP PARALLEL DEFAULT(NONE) &
! !$OMP PRIVATE(a, b, ab, i, j, ij, n) &
! !$OMP SHARED(nC, nOrb, nO, nS, pp_trip_Gam_B, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om)
! !$OMP DO COLLAPSE(2)
ab = 0
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
do n=1,nS
pp_Gam_C(ab,cd) = pp_Gam_C(ab,cd) &
- eh_rho(c,a,n)*eh_rho(b,d,n)/eh_Om(n) &
- eh_rho(a,c,n)*eh_rho(d,b,n)/eh_Om(n) &
+ eh_rho(d,a,n)*eh_rho(b,c,n)/eh_Om(n) &
+ eh_rho(a,d,n)*eh_rho(c,b,n)/eh_Om(n)
end do
end do
end do
end do
end do
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine
subroutine G_pp_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOO,nVV,eh_Om,eh_rho,pp_Gam_B)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOO,nVV
double precision,intent(in) :: eh_Om(nS)
double precision,intent(in) :: eh_rho(nOrb,nOrb,nS)
! 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)
! Initialization
pp_Gam_B(:,:) = 0d0
! !$OMP PARALLEL DEFAULT(NONE) &
! !$OMP PRIVATE(a, b, ab, i, j, ij, n) &
! !$OMP SHARED(nC, nOrb, nO, nS, pp_trip_Gam_B, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om)
! !$OMP DO COLLAPSE(2)
ab = 0
do a=nO+1,nOrb - nR
do b=a+1,nOrb - nR
ab = ab + 1
ij = 0
do i=nC+1,nO
do j=i+1,nO
ij = ij + 1
do n=1,nS
pp_Gam_B(ab,ij) = pp_Gam_B(ab,ij) &
- eh_rho(i,a,n)*eh_rho(b,j,n)/eh_Om(n) &
- eh_rho(a,i,n)*eh_rho(j,b,n)/eh_Om(n) &
+ eh_rho(j,a,n)*eh_rho(b,i,n)/eh_Om(n) &
+ eh_rho(a,j,n)*eh_rho(i,b,n)/eh_Om(n)
end do
end do
end do
end do
end do
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine

154
src/Parquet/G_screened_integrals.f90 Normal file
View File

@ -0,0 +1,154 @@
subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Gam,XpY,XmY,rho)
! Compute excitation densities
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nR,nS
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_Gam(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: XpY(nS,nS),XmY(nS,nS)
! Local variables
integer :: ia,jb,p,q,j,b
double precision :: X,Y
! Output variables
double precision,intent(out) :: rho(nOrb,nOrb,nS)
rho(:,:,:) = 0d0
! !$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
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))
rho(p,q,ia) = rho(p,q,ia) &
!+ (ERI(p,j,q,b) - ERI(p,j,b,q))*XpY(ia,jb) &
+ (ERI(p,j,q,b) - ERI(p,j,b,q))*X &
+ (ERI(p,b,q,j) - ERI(p,b,j,q))*Y &
+ 0d0*eh_Gam(p,j,q,b)*XpY(ia,jb)
end do
end do
end do
end do
end do
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine G_eh_screened_integral
subroutine G_pp_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_Gam,X1,Y1,rho1,X2,Y2,rho2)
! Compute excitation densities in the singlet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_Gam(nOrb,nOrb,nOrb,nOrb)
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: X1(nVV,nVV)
double precision,intent(in) :: Y1(nOO,nVV)
double precision,intent(in) :: X2(nVV,nOO)
double precision,intent(in) :: Y2(nOO,nOO)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
integer :: p,q
integer :: ab,cd,ij,kl
double precision,external :: Kronecker_delta
! Output variables
double precision,intent(out) :: rho1(nOrb,nOrb,nVV)
double precision,intent(out) :: rho2(nOrb,nOrb,nOO)
integer :: dim_1, dim_2
! Initialization
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)
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
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_Gam(p,q,c,d))*X1(cd,ab)
end do
end do
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho1(p,q,ab) = rho1(p,q,ab) &
+ (ERI(p,q,k,l) - ERI(p,q,l,k) + 1d0*pp_Gam(p,q,k,l))*Y1(kl,ab)
end do
end do
end do
end do
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
rho2(p,q,ij) = rho2(p,q,ij) &
+ (ERI(p,q,c,d) - ERI(p,q,d,c) + 1d0*pp_Gam(p,q,c,d))*X2(cd,ij)
end do
end do
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_Gam(p,q,k,l))*Y2(kl,ij)
end do
end do
end do
end do
end do
end do
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine G_pp_screened_integral

36
src/Parquet/RParquet.f90
View File

@ -9,8 +9,8 @@ subroutine RParquet(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
@ -262,14 +262,13 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
old_hh_sing_Om,hh_sing_rho,old_hh_trip_Om,hh_trip_rho, &
eh_sing_Gam_A)
call R_eh_singlet_Gamma_B(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_B)
if(.not.TDA) call R_eh_singlet_Gamma_B(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_B)
end if
end if
Aph(:,:) = Aph(:,:) + eh_sing_Gam_A(:,:)
Bph(:,:) = Bph(:,:) + eh_sing_Gam_B(:,:)
@ -320,11 +319,11 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
old_hh_sing_Om,hh_sing_rho,old_hh_trip_Om,hh_trip_rho, &
eh_trip_Gam_A)
call R_eh_triplet_Gamma_B(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_B)
if(.not.TDA) call R_eh_triplet_Gamma_B(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_B)
end if
@ -377,7 +376,7 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
else
call R_pp_singlet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,&
if(.not.TDA) call R_pp_singlet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,&
old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho,pp_sing_Gam_B)
call R_pp_singlet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVs,old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho,pp_sing_Gam_C)
call R_pp_singlet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOs,old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho,pp_sing_Gam_D)
@ -435,7 +434,7 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
else
call R_pp_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOt,nVVt,&
if(.not.TDA) call R_pp_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOt,nVVt,&
old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, pp_trip_Gam_B)
call R_pp_triplet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVt,old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, pp_trip_Gam_C)
call R_pp_triplet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOt,old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, pp_trip_Gam_D)
@ -503,7 +502,7 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
write(*,*) 'Computing singlet eh screened integrals...'
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 R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,sing_XpY,sing_XmY,eh_sing_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for singlet eh integrals =',t,' seconds'
@ -515,7 +514,7 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
write(*,*) 'Computing triplet eh screened integrals...'
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 R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,trip_XpY,trip_XmY,eh_trip_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for triplet eh integrals =',t,' seconds'
@ -549,7 +548,6 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
deallocate(X1t,Y1t,X2t,Y2t,pp_trip_Gam)
! Convergence criteria
err_2b = max(err_eh_sing,err_eh_trip,err_ee_sing,err_ee_trip,err_hh_sing,err_hh_trip)
call wall_time(end_2b)

19
src/Parquet/R_eh_singlet_Gam.f90
View File

@ -3,7 +3,6 @@ subroutine R_eh_singlet_Gamma(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
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
@ -40,8 +39,10 @@ subroutine R_eh_singlet_Gamma(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
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)
+ 0.5d0 * eh_sing_rho(s,p,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) &
+ 0.5d0 * eh_sing_rho(p,s,n)*eh_sing_rho(r,q,n)/eh_sing_Om(n) &
+ 1.5d0 * eh_trip_rho(s,p,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n) &
+ 1.5d0 * eh_trip_rho(p,s,n)*eh_trip_rho(r,q,n)/eh_trip_Om(n)
end do
do n=1,nVVs
@ -118,8 +119,10 @@ 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(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)
+ 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
@ -196,8 +199,10 @@ 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(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)
+ 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

26
src/Parquet/R_eh_triplet_Gam.f90
View File

@ -40,28 +40,30 @@ subroutine R_eh_triplet_Gamma(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
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)
+ 0.5d0 * eh_sing_rho(s,p,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) &
+ 0.5d0 * eh_sing_rho(p,s,n)*eh_sing_rho(r,q,n)/eh_sing_Om(n) &
- 0.5d0 * eh_trip_rho(s,p,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n) &
- 0.5d0 * eh_trip_rho(p,s,n)*eh_trip_rho(r,q,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)
- 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)
+ 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)
+ 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)
- hh_trip_rho(p,q,n) * hh_trip_rho(r,s,n)/hh_trip_Om(n)
end do
enddo
@ -118,8 +120,10 @@ 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(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)
+ 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
@ -196,8 +200,10 @@ 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(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)
+ 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

48
src/Parquet/R_pp_singlet_Gam.f90
View File

@ -33,10 +33,14 @@ 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(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)
- 0.5d0 * eh_sing_rho(r,p,n)*eh_sing_rho(q,s,n)/eh_sing_Om(n) &
- 0.5d0 * eh_sing_rho(p,r,n)*eh_sing_rho(s,q,n)/eh_sing_Om(n) &
+ 1.5d0 * eh_trip_rho(r,p,n)*eh_trip_rho(q,s,n)/eh_trip_Om(n) &
+ 1.5d0 * eh_trip_rho(p,r,n)*eh_trip_rho(s,q,n)/eh_trip_Om(n) &
- 0.5d0 * eh_sing_rho(s,p,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) &
- 0.5d0 * eh_sing_rho(p,s,n)*eh_sing_rho(r,q,n)/eh_sing_Om(n) &
+ 1.5d0 * eh_trip_rho(s,p,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n) &
+ 1.5d0 * eh_trip_rho(p,s,n)*eh_trip_rho(r,q,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 +95,14 @@ 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(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)
- 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) = pp_sing_Gam_D(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
@ -149,10 +157,14 @@ 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(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)
- 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)))
@ -207,10 +219,14 @@ 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(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)
- 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) = pp_sing_Gam_B(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))

52
src/Parquet/R_pp_triplet_Gam.f90
View File

@ -32,11 +32,16 @@ subroutine R_pp_triplet_Gamma(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh
do p = nC+1, nOrb-nR
do n=1,nS
pp_trip_Gam(p,q,r,s) = pp_trip_Gam(p,q,r,s) &
- 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)
- 0.5d0 * eh_sing_rho(r,p,n)*eh_sing_rho(q,s,n)/eh_sing_Om(n) &
- 0.5d0 * eh_sing_rho(p,r,n)*eh_sing_rho(s,q,n)/eh_sing_Om(n) &
- 0.5d0 * eh_trip_rho(r,p,n)*eh_trip_rho(q,s,n)/eh_trip_Om(n) &
- 0.5d0 * eh_trip_rho(p,r,n)*eh_trip_rho(s,q,n)/eh_trip_Om(n) &
+ 0.5d0 * eh_sing_rho(s,p,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) &
+ 0.5d0 * eh_sing_rho(p,s,n)*eh_sing_rho(r,q,n)/eh_sing_Om(n) &
+ 0.5d0 * eh_trip_rho(s,p,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n) &
+ 0.5d0 * eh_trip_rho(p,s,n)*eh_trip_rho(r,q,n)/eh_trip_Om(n)
end do
end do
@ -88,11 +93,16 @@ subroutine R_pp_triplet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOt,eh_sing_Om,eh_sing_rho,
kl = kl +1
do n=1,nS
pp_trip_Gam_D(ij,kl) = pp_trip_Gam_D(ij,kl) &
- 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)
- 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
end do
@ -144,11 +154,16 @@ subroutine R_pp_triplet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVt,eh_sing_Om,eh_sing_rho,
cd = cd +1
do n=1,nS
pp_trip_Gam_C(ab,cd) = pp_trip_Gam_C(ab,cd) &
- 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)
- 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
end do
@ -200,11 +215,16 @@ subroutine R_pp_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOt,nVVt,eh_sing_Om,eh_sing
ij = ij +1
do n=1,nS
pp_trip_Gam_B(ab,ij) = pp_trip_Gam_B(ab,ij) &
- 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)
- 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
end do

34
src/Parquet/R_screened_integrals.f90
View File

@ -1,4 +1,4 @@
subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,XpY,rho)
subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,XpY,XmY,rho)
! Compute excitation densities
implicit none
@ -7,10 +7,11 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,XpY,r
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)
double precision,intent(in) :: XpY(nS,nS),XmY(nS,nS)
! Local variables
integer :: ia,jb,p,q,j,b
double precision :: X,Y
! Output variables
double precision,intent(out) :: rho(nOrb,nOrb,nS)
@ -28,9 +29,14 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,XpY,r
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) &
+ 1d0*eh_sing_Gam(p,j,q,b)*XpY(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 &
+ 0d0*eh_sing_Gam(p,j,q,b)*XpY(ia,jb)
end do
end do
@ -42,7 +48,7 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,XpY,r
end subroutine R_eh_singlet_screened_integral
subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,XpY,rho)
subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,XpY,XmY,rho)
! Compute excitation densities
implicit none
@ -51,11 +57,12 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,XpY,r
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)
double precision,intent(in) :: XpY(nS,nS),XmY(nS,nS)
! Local variables
integer :: ia,jb,p,q,j,b
double precision :: X,Y
! Output variables
double precision,intent(out) :: rho(nOrb,nOrb,nS)
@ -72,9 +79,14 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,XpY,r
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) &
+ 1d0*eh_trip_Gam(p,j,q,b)*XpY(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 &
+ 0d0*eh_trip_Gam(p,j,q,b)*XpY(ia,jb)
end do
end do
end do

0
src/Parquet/amarie@HP-EliteBook-830.23187:1742801908 Normal file
View File

2
src/QuAcK/GQuAcK.f90
View File

@ -346,7 +346,7 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop
if(doParquet) then
call wall_time(start_Parquet)
call GParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body, &
nBas,nC,nO,nV,nR,nS, &
nBas2,nC,nO,nV,nR,nS, &
eHF,ERI_MO)
call wall_time(end_Parquet)