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

refactor in progress

Browse Source
This commit is contained in:
Antoine Marie 2025-03-25 17:20:02 +01:00
parent 1e0ed28f67
commit 9d8b33dbea
10 changed files with 649 additions and 1093 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

222
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 = .true.
logical :: print_ppLR = .true.
logical :: print_phLR = .false.
logical :: print_ppLR = .false.
! Input variables
@ -24,47 +24,50 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
integer :: n_it_1b,n_it_2b
double precision :: err_1b,err_2b
double precision :: err_eh, err_hh, err_ee
double precision :: err_eh, err_pp
double precision :: err_eig_eh, err_eig_hh, err_eig_ee
double precision :: start_t, end_t, t
double precision :: start_1b, end_1b, t_1b
double precision :: start_2b, end_2b, t_2b
integer :: nOO,nVV
! eh BSE
double precision :: EcRPA
double precision,allocatable :: Aph(:,:), Bph(:,:)
double precision,allocatable :: XpY(:,:),XmY(:,:)
double precision,allocatable :: XpY(:,:), XmY(:,:)
double precision,allocatable :: eh_Om(:), old_eh_Om(:)
double precision,allocatable :: eh_Gam_A(:,:),eh_Gam_B(:,:)
! pp BSE
double precision,allocatable :: Bpp(:,:), Cpp(:,:), Dpp(:,:)
double precision,allocatable :: X1(:,:),Y1(:,:)
double precision,allocatable :: ee_Om(:), old_ee_Om(:)
double precision,allocatable :: X2(:,:),Y2(:,:)
double precision,allocatable :: hh_Om(:), old_hh_Om(:)
double precision,allocatable :: eh_rho(:,:,:), ee_rho(:,:,:), hh_rho(:,:,:)
double precision,allocatable :: eh_Gam_A(:,:),eh_Gam_B(:,:)
double precision,allocatable :: pp_Gam_B(:,:),pp_Gam_C(:,:),pp_Gam_D(:,:)
double precision,allocatable :: eh_Gam(:,:,:,:),pp_Gam(:,:,:,:)
! Effective integrals
double precision,allocatable :: eh_rho(:,:,:), ee_rho(:,:,:), hh_rho(:,:,:)
! Reducible kernels
double precision,allocatable :: eh_Phi(:,:,:,:), pp_Phi(:,:,:,:)
double precision,allocatable :: old_eh_Phi(:,:,:,:), old_pp_Phi(:,:,:,:)
! One-body quantities
double precision,allocatable :: eQPlin(:),eQP(:),eOld(:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:)
double precision :: EcGM
! DIIS
integer :: max_diis,n_diis
double precision :: rcond
double precision,allocatable :: err_diis(:,:)
double precision,allocatable :: Om_diis(:,:)
double precision,allocatable :: err(:)
double precision,allocatable :: Om(:)
double precision :: alpha
! Output variables
! None
! Useful parameters
nOO = nO*(nO - 1)/2
nVV = nV*(nV - 1)/2
@ -72,15 +75,15 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
! DIIS parameters
max_diis = 10
n_diis = 0
rcond = 0d0
! max_diis = 10
! n_diis = 0
! rcond = 0d0
allocate(err_diis(nS+nOO+nVV,max_diis),Om_diis(nS+nOO+nVV,max_diis))
allocate(err(nS+nOO+nVV),Om(nS+nOO+nVV))
! allocate(err_diis(nS+nOO+nVV,max_diis),Om_diis(nS+nOO+nVV,max_diis))
! allocate(err(nS+nOO+nVV),Om(nS+nOO+nVV))
err_diis(:,:) = 0d0
Om_diis(:,:) = 0d0
! err_diis(:,:) = 0d0
! Om_diis(:,:) = 0d0
! Start
@ -115,6 +118,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
allocate(old_eh_Om(nS),old_ee_Om(nVV),old_hh_Om(nOO))
allocate(eh_rho(nOrb,nOrb,nS),ee_rho(nOrb,nOrb,nVV),hh_rho(nOrb,nOrb,nOO))
allocate(old_eh_Phi(nOrb,nOrb,nOrb,nOrb),old_pp_Phi(nOrb,nOrb,nOrb,nOrb))
! Initialization
@ -135,6 +139,9 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
old_ee_Om(:) = 0d0
old_hh_Om(:) = 0d0
old_eh_Phi(:,:,:,:) = 0d0
old_pp_Phi(:,:,:,:) = 0d0
!-----------------------------------------!
! Main loop for one-body self-consistency !
!-----------------------------------------!
@ -164,54 +171,6 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
write(*,*)' ***********************************'
write(*,*)
!--------------------------------!
! Compute effective interactions !
!--------------------------------!
! Memory allocation
allocate(eh_Gam(nOrb,nOrb,nOrb,nOrb))
allocate(pp_Gam(nOrb,nOrb,nOrb,nOrb))
! Build eh effective interaction
write(*,*) 'Computing eh effective interaction...'
if(n_it_2b == 1) then
eh_Gam(:,:,:,:) = 0d0
else
call wall_time(start_t)
call G_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_Gam)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A50,1X,F9.3,A8)') 'Wall time for eh Gamma =',t,' seconds'
write(*,*)
end if
! Build singlet pp effective interaction
write(*,*) 'Computing pp effective interaction...'
if(n_it_2b == 1) then
pp_Gam(:,:,:,:) = 0d0
else
call wall_time(start_t)
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(*,*)
end if
!-----------------!
! eh channel !
!-----------------!
@ -239,18 +198,14 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
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)
call G_eh_Gamma_A(nOrb,nC,nO,nR,nS,old_eh_Phi,old_pp_Phi,eh_Gam_A)
if(.not.TDA) call G_eh_Gamma_B(nOrb,nC,nO,nR,nS,old_eh_Phi,old_pp_Phi,eh_Gam_B)
end if
Aph(:,:) = Aph(:,:) + eh_Gam_A(:,:)
Bph(:,:) = Bph(:,:) + eh_Gam_B(:,:)
Aph(:,:) = Aph(:,:) + eh_Gam_A(:,:)
Bph(:,:) = Bph(:,:) + eh_Gam_B(:,:)
call phGLR(TDA,nS,Aph,Bph,EcRPA,eh_Om,XpY,XmY)
@ -262,7 +217,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
if(print_phLR) call print_excitation_energies('phBSE@Parquet','eh generalized',nS,eh_Om)
err_eh = maxval(abs(old_eh_Om - eh_Om))
err_eig_eh = maxval(abs(old_eh_Om - eh_Om))
deallocate(Aph,Bph,eh_Gam_A,eh_Gam_B)
@ -297,9 +252,9 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
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)
if(.not.TDA) call G_pp_Gamma_B(nOrb,nC,nO,nR,nOO,nVV,old_eh_Phi,pp_Gam_B)
call G_pp_Gamma_C(nOrb,nO,nR,nVV,old_eh_Phi,pp_Gam_C)
call G_pp_Gamma_D(nOrb,nC,nO,nOO,old_eh_Phi,pp_Gam_D)
end if
@ -317,17 +272,17 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
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))
err_eig_ee = maxval(abs(old_ee_Om - ee_Om))
err_eig_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(*,*) ' Two-body (eigenvalue) 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(*,'(1X,A30,F10.6)')'Error for eh channel = ',err_eig_eh
write(*,'(1X,A30,F10.6)')'Error for pp channel = ',max(err_eig_ee,err_eig_hh)
write(*,*) '----------------------------------------'
write(*,*)
@ -335,24 +290,24 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
! DIIS extrapolation !
!--------------------!
err( 1:nS ) = eh_Om(:) - old_eh_Om(:)
err( nS+1:nS+nVV ) = ee_Om(:) - old_ee_Om(:)
err(nVV+nS+1:nS+nVV+nOO) = hh_Om(:) - old_hh_Om(:)
! err( 1:nS ) = eh_Om(:) - old_eh_Om(:)
! err( nS+1:nS+nVV ) = ee_Om(:) - old_ee_Om(:)
! err(nVV+nS+1:nS+nVV+nOO) = hh_Om(:) - old_hh_Om(:)
Om( 1:nS ) = eh_Om(:)
Om( nS+1:nS+nVV ) = ee_Om(:)
Om(nVV+nS+1:nS+nVV+nOO) = hh_Om(:)
! Om( 1:nS ) = eh_Om(:)
! Om( nS+1:nS+nVV ) = ee_Om(:)
! Om(nVV+nS+1:nS+nVV+nOO) = hh_Om(:)
if(max_diis > 1) then
! if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nS+nOO+nVV,nS+nOO+nVV,n_diis,err_diis,Om_diis,err,Om)
! n_diis = min(n_diis+1,max_diis)
! call DIIS_extrapolation(rcond,nS+nOO+nVV,nS+nOO+nVV,n_diis,err_diis,Om_diis,err,Om)
end if
! end if
eh_Om(:) = Om( 1:nS )
ee_Om(:) = Om( nS+1:nS+nVV )
hh_Om(:) = Om(nVV+nS+1:nS+nVV+nOO)
! eh_Om(:) = Om( 1:nS )
! ee_Om(:) = Om( nS+1:nS+nVV )
! hh_Om(:) = Om(nVV+nS+1:nS+nVV+nOO)
!----------!
! Updating !
@ -379,35 +334,90 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
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 G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,old_eh_Phi,old_pp_Phi,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)
deallocate(XpY,XmY)
! 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 G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,old_eh_Phi,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)
deallocate(X1,Y1,X2,Y2)
!----------------------------!
! Compute reducible kernels !
!----------------------------!
! Memory allocation
allocate(eh_Phi(nOrb,nOrb,nOrb,nOrb))
allocate(pp_Phi(nOrb,nOrb,nOrb,nOrb))
! Build eh reducible kernels
write(*,*) 'Computing eh reducible kernel ...'
call wall_time(start_t)
call G_eh_Phi(nOrb,nC,nR,nS,old_eh_Om,eh_rho,eh_Phi)
call wall_time(end_t)
t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for eh reducible kernel =',t,' seconds'
write(*,*)
! Build pp reducible kernels
write(*,*) 'Computing pp reducible kernel ...'
call wall_time(start_t)
call G_pp_Phi(nOrb,nC,nR,nOO,nVV,old_ee_Om,ee_rho,old_hh_Om,hh_rho,pp_Phi)
call wall_time(end_t)
t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for pp reducible kernel =',t,' seconds'
write(*,*)
! alpha = 0.01d0
! eh_Phi(:,:,:,:) = alpha * eh_Phi(:,:,:,:) + (1d0 - alpha) * old_eh_Phi(:,:,:,:)
! pp_Phi(:,:,:,:) = alpha * pp_Phi(:,:,:,:) + (1d0 - alpha) * old_pp_Phi(:,:,:,:)
err_eh = maxval(abs(old_eh_Phi - eh_Phi))
err_pp = maxval(abs(old_pp_Phi - pp_Phi))
old_eh_Phi(:,:,:,:) = eh_Phi(:,:,:,:)
old_pp_Phi(:,:,:,:) = pp_Phi(:,:,:,:)
! Free memory
deallocate(eh_Phi,pp_Phi)
!--------------------!
! DIIS extrapolation !
!--------------------!
write(*,*) '----------------------------------------'
write(*,*) ' Two-body (kernel) convergence '
write(*,*) '----------------------------------------'
write(*,'(1X,A30,F10.6)')'Error for eh channel = ',err_eh
write(*,'(1X,A30,F10.6)')'Error for pp channel = ',err_pp
write(*,*) '----------------------------------------'
write(*,*)
! Convergence criteria
err_2b = max(err_eh,err_ee,err_hh)
err_2b = max(err_eh,err_pp)
call wall_time(end_2b)
t_2b = end_2b - start_2b
write(*,'(A50,1X,F9.3,A8)') 'Wall time for two-body iteration =',t_2b,' seconds'
write(*,*)
end do
!---------------------------------------------!
! End main loop for two-body self-consistency !

121
src/Parquet/G_eh_Gam.f90
View File

@ -1,79 +1,16 @@
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)
subroutine G_eh_Gamma_A(nOrb,nC,nO,nR,nS,eh_Phi,pp_Phi,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 :: 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_Gam(p,q,r,s) = eh_Gam(p,q,r,s) &
+ 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
do n=1,nVV
eh_Gam(p,q,r,s) = eh_Gam(p,q,r,s) &
+ 2d0 * ee_rho(p,q,n)*ee_rho(r,s,n)/ee_Om(n)
end do
do n=1,nOO
eh_Gam(p,q,r,s) = eh_Gam(p,q,r,s) &
- 2d0 * hh_rho(p,q,n)*hh_rho(r,s,n)/hh_Om(n)
end do
enddo
enddo
enddo
enddo
end subroutine G_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)
integer,intent(in) :: nOrb,nC,nO,nR,nS
double precision,intent(in) :: eh_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_Phi(nOrb,nOrb,nOrb,nOrb)
! 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)
@ -91,21 +28,7 @@ subroutine G_eh_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOO,nVV, &
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
eh_Gam_A(ia,jb) = - eh_Phi(a,j,b,i) + pp_Phi(a,j,i,b)
enddo
enddo
@ -114,27 +37,19 @@ subroutine G_eh_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOO,nVV, &
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)
subroutine G_eh_Gamma_B(nOrb,nC,nO,nR,nS,eh_Phi,pp_Phi,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)
integer,intent(in) :: nOrb,nC,nO,nR,nS
double precision,intent(in) :: eh_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_Phi(nOrb,nOrb,nOrb,nOrb)
! 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)
@ -152,21 +67,7 @@ subroutine G_eh_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOO,nVV, &
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
eh_Gam_B(ia,jb) = - eh_Phi(a,b,j,i) + pp_Phi(a,b,i,j)
enddo
enddo

103
src/Parquet/G_pp_Gam.f90
View File

@ -1,66 +1,15 @@
subroutine G_pp_Gamma(nOrb,nC,nO,nV,nR,nS,eh_Om,eh_rho,pp_Gam)
subroutine G_pp_Gamma_D(nOrb,nC,nO,nOO,eh_Phi,pp_Gam_D)
! 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)
integer,intent(in) :: nOrb,nC,nO,nOO
double precision,intent(in) :: eh_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: i,j,k,l
integer :: ij,kl
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_Gam_D(nOO,nOO)
@ -82,15 +31,9 @@ subroutine G_pp_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOO,eh_Om,eh_rho,pp_Gam_D)
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
pp_Gam_D(ij,kl) = eh_Phi(i,j,k,l) - eh_Phi(i,j,l,k)
end do
end do
end do
@ -100,21 +43,18 @@ subroutine G_pp_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOO,eh_Om,eh_rho,pp_Gam_D)
end subroutine
subroutine G_pp_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVV,eh_Om,eh_rho,pp_Gam_C)
subroutine G_pp_Gamma_C(nOrb,nO,nR,nVV,eh_Phi,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)
integer,intent(in) :: nOrb,nO,nR,nVV
double precision,intent(in) :: eh_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: a,b,c,d
integer :: ab,cd
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_Gam_C(nVV,nVV)
@ -137,15 +77,7 @@ subroutine G_pp_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVV,eh_Om,eh_rho,pp_Gam_C)
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
pp_Gam_C(ab,cd) = eh_Phi(a,b,c,d) - eh_Phi(a,b,d,c)
end do
end do
@ -156,15 +88,14 @@ subroutine G_pp_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVV,eh_Om,eh_rho,pp_Gam_C)
end subroutine
subroutine G_pp_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOO,nVV,eh_Om,eh_rho,pp_Gam_B)
subroutine G_pp_Gamma_B(nOrb,nC,nO,nR,nOO,nVV,eh_Phi,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)
integer,intent(in) :: nOrb,nC,nO,nR,nOO,nVV
double precision,intent(in) :: eh_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: a,b,i,j
@ -193,13 +124,7 @@ subroutine G_pp_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOO,nVV,eh_Om,eh_rho,pp_Gam_B)
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
pp_Gam_B(ab,ij) = eh_Phi(a,b,i,j) - eh_Phi(a,b,j,i)
end do
end do

38
src/Parquet/G_screened_integrals.f90
View File

@ -1,4 +1,4 @@
subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Gam,XpY,XmY,rho)
subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho)
! Compute excitation densities
implicit none
@ -6,7 +6,8 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Gam,XpY,XmY,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_Gam(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: XpY(nS,nS),XmY(nS,nS)
! Local variables
@ -35,13 +36,10 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Gam,XpY,XmY,rho)
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 &
+ 1d0*eh_Gam(p,j,q,b)*X &
+ 1d0*eh_Gam(p,b,q,j)*Y
rho(p,q,ia) = (ERI(q,j,p,b) - ERI(q,j,b,p)) * X &
+ (- eh_Phi(q,j,b,p) + pp_Phi(q,j,p,b)) * X &
+ (ERI(q,b,p,j) - ERI(q,b,j,p)) * Y &
+ (- eh_Phi(q,b,j,p) + pp_Phi(q,b,p,j)) * Y
end do
end do
@ -53,7 +51,7 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Gam,XpY,XmY,rho)
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)
subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2,Y2,rho2)
! Compute excitation densities in the singlet pp channel
@ -61,9 +59,9 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_Gam,X1,Y1,rho1
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR
integer,intent(in) :: nOrb,nC,nO,nR
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_Gam(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_Phi(nOrb,nOrb,nOrb,nOrb)
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: X1(nVV,nVV)
@ -107,8 +105,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_Gam,X1,Y1,rho1
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)
rho1(p,q,ab) = ( ERI(p,q,c,d) - ERI(p,q,d,c) ) * X1(cd,ab) &
+ ( eh_Phi(p,q,c,d) - eh_Phi(p,q,d,c) ) * X1(cd,ab)
end do
end do
@ -116,8 +114,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_Gam,X1,Y1,rho1
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)
rho1(p,q,ab) = ( ERI(p,q,k,l) - ERI(p,q,l,k) ) * Y1(kl,ab) &
+ ( eh_Phi(p,q,k,l) - eh_Phi(p,q,l,k) ) * Y1(kl,ab)
end do
end do
@ -132,8 +130,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_Gam,X1,Y1,rho1
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)
rho2(p,q,ij) = ( ERI(p,q,c,d) - ERI(p,q,d,c) ) * X2(cd,ij) &
+ ( eh_Phi(p,q,c,d) - eh_Phi(p,q,d,c) ) * X2(cd,ij)
end do
end do
@ -141,8 +139,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_Gam,X1,Y1,rho1
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)
rho2(p,q,ij) = ( ERI(p,q,k,l) - ERI(p,q,l,k) ) * Y2(kl,ij) &
+ ( eh_Phi(p,q,k,l) - eh_Phi(p,q,l,k) ) * Y2(kl,ij)
end do
end do
end do

271
src/Parquet/RParquet.f90
View File

@ -26,23 +26,28 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
integer :: n_it_1b,n_it_2b
double precision :: err_1b,err_2b
double precision :: err_eh_sing,err_eh_trip
double precision :: err_hh_sing,err_hh_trip
double precision :: err_ee_sing,err_ee_trip
double precision :: err_eig_eh_sing,err_eig_eh_trip
double precision :: err_eig_hh_sing,err_eig_hh_trip
double precision :: err_eig_ee_sing,err_eig_ee_trip
double precision :: err_eh_sing, err_eh_trip
double precision :: err_pp_sing, err_pp_trip
double precision :: start_t, end_t, t
double precision :: start_1b, end_1b, t_1b
double precision :: start_2b, end_2b, t_2b
integer :: nOOs,nOOt
integer :: nVVs,nVVt
! eh BSE
double precision :: EcRPA
double precision,allocatable :: Aph(:,:), Bph(:,:)
double precision,allocatable :: sing_XpY(:,:),trip_XpY(:,:)
double precision,allocatable :: sing_XmY(:,:),trip_XmY(:,:)
double precision,allocatable :: eh_sing_Om(:), old_eh_sing_Om(:)
double precision,allocatable :: eh_trip_Om(:), old_eh_trip_Om(:)
double precision,allocatable :: eh_sing_Gam_A(:,:),eh_sing_Gam_B(:,:)
double precision,allocatable :: eh_trip_Gam_A(:,:),eh_trip_Gam_B(:,:)
! pp BSE
double precision,allocatable :: Bpp(:,:), Cpp(:,:), Dpp(:,:)
double precision,allocatable :: X1s(:,:),X1t(:,:)
double precision,allocatable :: Y1s(:,:),Y1t(:,:)
@ -52,18 +57,18 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
double precision,allocatable :: Y2s(:,:),Y2t(:,:)
double precision,allocatable :: hh_sing_Om(:), old_hh_sing_Om(:)
double precision,allocatable :: hh_trip_Om(:), old_hh_trip_Om(:)
double precision,allocatable :: pp_sing_Gam_B(:,:),pp_sing_Gam_C(:,:),pp_sing_Gam_D(:,:)
double precision,allocatable :: pp_trip_Gam_B(:,:),pp_trip_Gam_C(:,:),pp_trip_Gam_D(:,:)
! Effective integrals
double precision,allocatable :: eh_sing_rho(:,:,:),eh_trip_rho(:,:,:)
double precision,allocatable :: ee_sing_rho(:,:,:),hh_sing_rho(:,:,:)
double precision,allocatable :: ee_trip_rho(:,:,:),hh_trip_rho(:,:,:)
double precision,allocatable :: eh_sing_Gam_A(:,:),eh_sing_Gam_B(:,:)
double precision,allocatable :: eh_trip_Gam_A(:,:),eh_trip_Gam_B(:,:)
double precision,allocatable :: pp_sing_Gam_B(:,:),pp_sing_Gam_C(:,:),pp_sing_Gam_D(:,:)
double precision,allocatable :: pp_trip_Gam_B(:,:),pp_trip_Gam_C(:,:),pp_trip_Gam_D(:,:)
double precision,allocatable :: eh_sing_Gam(:,:,:,:),eh_trip_Gam(:,:,:,:)
double precision,allocatable :: pp_sing_Gam(:,:,:,:),pp_trip_Gam(:,:,:,:)
! Reducible kernels
double precision,allocatable :: eh_sing_Phi(:,:,:,:), eh_trip_Phi(:,:,:,:)
double precision,allocatable :: old_eh_sing_Phi(:,:,:,:), old_eh_trip_Phi(:,:,:,:)
double precision,allocatable :: pp_sing_Phi(:,:,:,:), pp_trip_Phi(:,:,:,:)
double precision,allocatable :: old_pp_sing_Phi(:,:,:,:), old_pp_trip_Phi(:,:,:,:)
! One-body quantities
double precision,allocatable :: eQPlin(:),eQP(:),eOld(:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:)
@ -71,7 +76,8 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
! Output variables
! None
! Useful parameters
nOOs = nO*(nO + 1)/2
nVVs = nV*(nV + 1)/2
nOOt = nO*(nO - 1)/2
@ -114,6 +120,8 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
allocate(eh_sing_rho(nOrb,nOrb,nS),eh_trip_rho(nOrb,nOrb,nS))
allocate(ee_sing_rho(nOrb,nOrb,nVVs),hh_sing_rho(nOrb,nOrb,nOOs))
allocate(ee_trip_rho(nOrb,nOrb,nVVt),hh_trip_rho(nOrb,nOrb,nOOt))
allocate(old_eh_sing_Phi(nOrb,nOrb,nOrb,nOrb),old_eh_trip_Phi(nOrb,nOrb,nOrb,nOrb))
allocate(old_pp_sing_Phi(nOrb,nOrb,nOrb,nOrb),old_pp_trip_Phi(nOrb,nOrb,nOrb,nOrb))
! Initialization
@ -139,6 +147,11 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
old_ee_trip_Om(:) = 1d0
old_hh_sing_Om(:) = 1d0
old_hh_trip_Om(:) = 1d0
old_eh_sing_Phi(:,:,:,:) = 0d0
old_eh_trip_Phi(:,:,:,:) = 0d0
old_pp_sing_Phi(:,:,:,:) = 0d0
old_pp_trip_Phi(:,:,:,:) = 0d0
!-----------------------------------------!
! Main loop for one-body self-consistency !
@ -169,66 +182,6 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
write(*,*)' ***********************************'
write(*,*)
!--------------------------------!
! Compute effective interactions !
!--------------------------------!
! Memory allocation
allocate(eh_sing_Gam(nOrb,nOrb,nOrb,nOrb))
allocate(eh_trip_Gam(nOrb,nOrb,nOrb,nOrb))
allocate(pp_sing_Gam(nOrb,nOrb,nOrb,nOrb))
allocate(pp_trip_Gam(nOrb,nOrb,nOrb,nOrb))
! Build singlet eh effective interaction
write(*,*) 'Computing singlet eh effective interaction...'
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(*,'(A50,1X,F9.3,A8)') 'Wall time for eh singlet Gamma =',t,' seconds'
write(*,*)
! Build triplet eh effective interaction
write(*,*) 'Computing triplet eh effective interaction...'
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(*,'(A50,1X,F9.3,A8)') 'Wall time for eh triplet Gamma =',t,' seconds'
write(*,*)
! Build singlet pp effective interaction
write(*,*) 'Computing singlet pp effective interaction...'
call wall_time(start_t)
call R_pp_singlet_Gamma(nOrb,nC,nR,nS,old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho,pp_sing_Gam)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A50,1X,F9.3,A8)') 'Wall time for pp singlet Gamma =',t,' seconds'
write(*,*)
! Build triplet pp effective interaction
write(*,*) 'Computing triplet pp effective interaction...'
call wall_time(start_t)
call R_pp_triplet_Gamma(nOrb,nC,nR,nS,old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho,pp_trip_Gam)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A50,1X,F9.3,A8)') 'Wall time for pp triplet Gamma =',t,' seconds'
write(*,*)
!-----------------!
! Density channel !
!-----------------!
@ -256,17 +209,13 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
else
call R_eh_singlet_Gamma_A(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_A)
call R_eh_singlet_Gamma_A(nOrb,nC,nO,nR,nS, &
old_eh_sing_Phi,old_eh_trip_Phi,old_pp_sing_Phi,old_pp_trip_Phi, &
eh_sing_Gam_A)
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)
if(.not.TDA) call R_eh_singlet_Gamma_B(nOrb,nC,nO,nR,nS, &
old_eh_sing_Phi,old_eh_trip_Phi,old_pp_sing_Phi,old_pp_trip_Phi, &
eh_sing_Gam_B)
end if
Aph(:,:) = Aph(:,:) + eh_sing_Gam_A(:,:)
@ -282,7 +231,7 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
if(print_phLR) call print_excitation_energies('phBSE@Parquet','singlet',nS,eh_sing_Om)
err_eh_sing = maxval(abs(old_eh_sing_Om - eh_sing_Om))
err_eig_eh_sing = maxval(abs(old_eh_sing_Om - eh_sing_Om))
deallocate(Aph,Bph,eh_sing_Gam_A,eh_sing_Gam_B)
@ -313,17 +262,13 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
else
call R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,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_A)
call R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS, &
old_eh_sing_Phi,old_eh_trip_Phi,old_pp_sing_Phi,old_pp_trip_Phi, &
eh_trip_Gam_A)
if(.not.TDA) call R_eh_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,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, &
old_eh_sing_Phi,old_eh_trip_Phi,old_pp_sing_Phi,old_pp_trip_Phi, &
eh_trip_Gam_B)
end if
@ -340,7 +285,7 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
if(print_phLR) call print_excitation_energies('phBSE@Parquet','triplet',nS,eh_trip_Om)
err_eh_trip = maxval(abs(old_eh_trip_Om - eh_trip_Om))
err_eig_eh_trip = maxval(abs(old_eh_trip_Om - eh_trip_Om))
deallocate(Aph,Bph,eh_trip_Gam_A,eh_trip_Gam_B)
@ -376,10 +321,10 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
else
if(.not.TDA) call R_pp_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)
if(.not.TDA) call R_pp_singlet_Gamma_B(nOrb,nC,nO,nR,nOOs,nVVs,&
old_eh_sing_Phi,old_eh_trip_Phi,pp_sing_Gam_B)
call R_pp_singlet_Gamma_C(nOrb,nO,nR,nVVs,old_eh_sing_Phi,old_eh_trip_Phi,pp_sing_Gam_C)
call R_pp_singlet_Gamma_D(nOrb,nC,nO,nOOs,old_eh_sing_Phi,old_eh_trip_Phi,pp_sing_Gam_D)
end if
@ -397,8 +342,8 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
if(print_ppLR) call print_excitation_energies('ppBSE@Parquet','2p (singlets)',nVVs,ee_sing_Om)
if(print_ppLR) call print_excitation_energies('ppBSE@Parquet','2h (singlets)',nOOs,hh_sing_Om)
err_ee_sing = maxval(abs(old_ee_sing_Om - ee_sing_Om))
err_hh_sing = maxval(abs(old_hh_sing_Om - hh_sing_Om))
err_eig_ee_sing = maxval(abs(old_ee_sing_Om - ee_sing_Om))
err_eig_hh_sing = maxval(abs(old_hh_sing_Om - hh_sing_Om))
deallocate(Bpp,Cpp,Dpp,pp_sing_Gam_B,pp_sing_Gam_C,pp_sing_Gam_D)
@ -434,10 +379,10 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
else
if(.not.TDA) call R_pp_triplet_Gamma_B(nOrb,nC,nO,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)
if(.not.TDA) call R_pp_triplet_Gamma_B(nOrb,nC,nO,nR,nS,nOOt,nVVt,&
old_eh_sing_Phi,old_eh_trip_Phi,pp_trip_Gam_B)
call R_pp_triplet_Gamma_C(nOrb,nO,nR,nS,nVVt,old_eh_sing_Phi,old_eh_trip_Phi,pp_trip_Gam_C)
call R_pp_triplet_Gamma_D(nOrb,nC,nO,nS,nOOt,old_eh_sing_Phi,old_eh_trip_Phi,pp_trip_Gam_D)
end if
@ -456,18 +401,18 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
if(print_ppLR) call print_excitation_energies('ppBSE@Parquet','2p (triplets)',nVVt,ee_trip_Om)
if(print_ppLR) call print_excitation_energies('ppBSE@Parquet','2h (triplets)',nOOt,hh_trip_Om)
err_ee_trip = maxval(abs(old_ee_trip_Om - ee_trip_Om))
err_hh_trip = maxval(abs(old_hh_trip_Om - hh_trip_Om))
err_eig_ee_trip = maxval(abs(old_ee_trip_Om - ee_trip_Om))
err_eig_hh_trip = maxval(abs(old_hh_trip_Om - hh_trip_Om))
deallocate(Bpp,Cpp,Dpp,pp_trip_Gam_B,pp_trip_Gam_C,pp_trip_Gam_D)
write(*,*) '----------------------------------------'
write(*,*) ' Two-body convergence '
write(*,*) '----------------------------------------'
write(*,'(1X,A30,F10.6)')'Error for density channel = ',err_eh_sing
write(*,'(1X,A30,F10.6)')'Error for magnetic channel = ',err_eh_trip
write(*,'(1X,A30,F10.6)')'Error for singlet channel = ',max(err_ee_sing,err_hh_sing)
write(*,'(1X,A30,F10.6)')'Error for triplet channel = ',max(err_ee_trip,err_hh_trip)
write(*,'(1X,A30,F10.6)')'Error for density channel = ',err_eig_eh_sing
write(*,'(1X,A30,F10.6)')'Error for magnetic channel = ',err_eig_eh_trip
write(*,'(1X,A30,F10.6)')'Error for singlet channel = ',max(err_eig_ee_sing,err_eig_hh_sing)
write(*,'(1X,A30,F10.6)')'Error for triplet channel = ',max(err_eig_ee_trip,err_eig_hh_trip)
write(*,*) '----------------------------------------'
write(*,*)
@ -502,53 +447,135 @@ 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,sing_XmY,eh_sing_rho)
call R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,old_eh_sing_Phi,old_eh_trip_Phi,old_pp_sing_Phi,old_pp_trip_Phi, &
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'
write(*,*)
! Done with eigenvectors and kernel
deallocate(sing_XpY,sing_XmY,eh_sing_Gam)
deallocate(sing_XpY,sing_XmY)
! Build triplet eh screened integrals
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,trip_XmY,eh_trip_rho)
call R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,old_eh_sing_Phi,old_eh_trip_Phi,old_pp_sing_Phi,old_pp_trip_Phi, &
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'
write(*,*)
! Done with eigenvectors and kernel
deallocate(trip_XpY,trip_XmY,eh_trip_Gam)
deallocate(trip_XpY,trip_XmY)
! Build singlet pp screened integrals
write(*,*) 'Computing singlet pp screened integrals...'
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 R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOOs,nVVs,ERI,old_eh_sing_Phi,old_eh_trip_Phi, &
X1s,Y1s,ee_sing_rho,X2s,Y2s,hh_sing_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 singlet pp integrals =',t,' seconds'
write(*,*)
deallocate(X1s,Y1s,X2s,Y2s,pp_sing_Gam)
deallocate(X1s,Y1s,X2s,Y2s)
! Build triplet pp screened integrals
write(*,*) 'Computing triplet pp screened integrals...'
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 R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOOt,nVVt,ERI,old_eh_sing_Phi,old_eh_trip_Phi, &
X1t,Y1t,ee_trip_rho,X2t,Y2t,hh_trip_rho)
call wall_time(end_t)
t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for triplet pp integrals =',t,' seconds'
write(*,*)
! Done with eigenvectors and kernel
deallocate(X1t,Y1t,X2t,Y2t,pp_trip_Gam)
deallocate(X1t,Y1t,X2t,Y2t)
!----------------------------!
! Compute reducible kernels !
!----------------------------!
! Memory allocation
allocate(eh_sing_Phi(nOrb,nOrb,nOrb,nOrb))
allocate(eh_trip_Phi(nOrb,nOrb,nOrb,nOrb))
allocate(pp_sing_Phi(nOrb,nOrb,nOrb,nOrb))
allocate(pp_trip_Phi(nOrb,nOrb,nOrb,nOrb))
! Build singlet eh reducible kernels
write(*,*) 'Computing singlet eh reducible kernel ...'
call wall_time(start_t)
call R_eh_singlet_Phi(nOrb,nC,nR,nS,old_eh_sing_Om,eh_sing_rho,eh_sing_Phi)
call wall_time(end_t)
t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for singlet eh reducible kernel =',t,' seconds'
write(*,*)
! Build triplet eh reducible kernels
write(*,*) 'Computing triplet eh reducible kernel ...'
call wall_time(start_t)
call R_eh_triplet_Phi(nOrb,nC,nR,nS,old_eh_trip_Om,eh_trip_rho,eh_trip_Phi)
call wall_time(end_t)
t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for triplet eh reducible kernel =',t,' seconds'
write(*,*)
! Build singlet pp reducible kernels
write(*,*) 'Computing singlet pp reducible kernel ...'
call wall_time(start_t)
call R_pp_singlet_Phi(nOrb,nC,nR,nOOs,nVVs,old_ee_sing_Om,ee_sing_rho,old_hh_sing_Om,hh_sing_rho,pp_sing_Phi)
call wall_time(end_t)
t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for singlet pp reducible kernel =',t,' seconds'
write(*,*)
! Build triplet pp reducible kernels
write(*,*) 'Computing triplet pp reducible kernel ...'
call wall_time(start_t)
call R_pp_triplet_Phi(nOrb,nC,nR,nOOt,nVVt,old_ee_trip_Om,ee_trip_rho,old_hh_trip_Om,hh_trip_rho,pp_trip_Phi)
call wall_time(end_t)
t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for triplet pp reducible kernel =',t,' seconds'
write(*,*)
err_eh_sing = maxval(abs(old_eh_sing_Phi - eh_sing_Phi))
err_eh_trip = maxval(abs(old_eh_trip_Phi - eh_trip_Phi))
err_pp_sing = maxval(abs(old_pp_sing_Phi - pp_sing_Phi))
err_pp_trip = maxval(abs(old_pp_trip_Phi - pp_trip_Phi))
old_eh_sing_Phi(:,:,:,:) = eh_sing_Phi(:,:,:,:)
old_eh_trip_Phi(:,:,:,:) = eh_trip_Phi(:,:,:,:)
old_pp_sing_Phi(:,:,:,:) = pp_sing_Phi(:,:,:,:)
old_pp_trip_Phi(:,:,:,:) = pp_trip_Phi(:,:,:,:)
! Free memory
deallocate(eh_sing_Phi,eh_trip_Phi,pp_sing_Phi,pp_trip_Phi)
!--------------------!
! DIIS extrapolation !
!--------------------!
write(*,*) '----------------------------------------'
write(*,*) ' Two-body (kernel) convergence '
write(*,*) '----------------------------------------'
write(*,'(1X,A30,F10.6)')'Error for singlet eh channel = ',err_eh_sing
write(*,'(1X,A30,F10.6)')'Error for triplet eh channel = ',err_eh_trip
write(*,'(1X,A30,F10.6)')'Error for singlet pp channel = ',err_pp_sing
write(*,'(1X,A30,F10.6)')'Error for triplet pp channel = ',err_pp_trip
write(*,*) '----------------------------------------'
write(*,*)
! Convergence criteria
err_2b = max(err_eh_sing,err_eh_trip,err_ee_sing,err_ee_trip,err_hh_sing,err_hh_trip)
err_2b = max(err_eh_sing,err_eh_trip,err_pp_sing,err_pp_trip)
call wall_time(end_2b)
t_2b = end_2b - start_2b

214
src/Parquet/R_eh_singlet_Gam.f90
View File

@ -1,106 +1,20 @@
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) &
+ 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
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, &
hh_sing_Om,hh_sing_rho,hh_trip_Om,hh_trip_rho, eh_sing_Gam_A)
subroutine R_eh_singlet_Gamma_A(nOrb,nC,nO,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing_Phi,pp_trip_Phi,eh_sing_Gam_A)
! 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)
integer,intent(in) :: nOrb,nC,nO,nR,nS
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_trip_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: i,a,j,b
integer :: ia,jb
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: eh_sing_Gam_A(nS,nS)
@ -117,33 +31,33 @@ subroutine R_eh_singlet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
do b=nO+1,norb-nR
jb = jb + 1
do n=1,nS
eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) &
+ 0.5d0 * eh_sing_rho(b,a,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) &
+ 0.5d0 * eh_sing_rho(a,b,n)*eh_sing_rho(i,j,n)/eh_sing_Om(n) &
+ 1.5d0 * eh_trip_rho(b,a,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n) &
+ 1.5d0 * eh_trip_rho(a,b,n)*eh_trip_rho(i,j,n)/eh_trip_Om(n)
end do
! do n=1,nS
! eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) &
! + 0.5d0 * eh_sing_rho(b,a,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) &
! + 0.5d0 * eh_sing_rho(a,b,n)*eh_sing_rho(i,j,n)/eh_sing_Om(n) &
! + 1.5d0 * eh_trip_rho(b,a,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n) &
! + 1.5d0 * eh_trip_rho(a,b,n)*eh_trip_rho(i,j,n)/eh_trip_Om(n)
! end do
do n=1,nVVs
eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) &
+ ee_sing_rho(a,j,n)*ee_sing_rho(i,b,n)/ee_sing_Om(n)
end do
! do n=1,nVVs
! eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) &
! + ee_sing_rho(a,j,n)*ee_sing_rho(i,b,n)/ee_sing_Om(n)
! end do
do n=1,nOOs
eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) &
- hh_sing_rho(a,j,n)*hh_sing_rho(i,b,n)/hh_sing_Om(n)
end do
! do n=1,nOOs
! eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) &
! - hh_sing_rho(a,j,n)*hh_sing_rho(i,b,n)/hh_sing_Om(n)
! end do
do n=1,nVVt
eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) &
+ 3d0 * ee_trip_rho(a,j,n)*ee_trip_rho(i,b,n)/ee_trip_Om(n)
end do
! do n=1,nVVt
! eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) &
! + 3d0 * ee_trip_rho(a,j,n)*ee_trip_rho(i,b,n)/ee_trip_Om(n)
! end do
do n=1,nOOt
eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) &
- 3d0 * hh_trip_rho(a,j,n)*hh_trip_rho(i,b,n)/hh_trip_Om(n)
end do
! do n=1,nOOt
! eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) &
! - 3d0 * hh_trip_rho(a,j,n)*hh_trip_rho(i,b,n)/hh_trip_Om(n)
! end do
enddo
enddo
@ -152,35 +66,23 @@ subroutine R_eh_singlet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
end subroutine R_eh_singlet_Gamma_A
subroutine R_eh_singlet_Gamma_B(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_B)
subroutine R_eh_singlet_Gamma_B(nOrb,nC,nO,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing_Phi,pp_trip_Phi,eh_sing_Gam_B)
! 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)
integer,intent(in) :: nOrb,nC,nO,nR,nS
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_trip_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: i,a,j,b
integer :: ia,jb
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: eh_sing_Gam_B(nS,nS)
@ -197,33 +99,33 @@ subroutine R_eh_singlet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
do b=nO+1,norb-nR
jb = jb + 1
do n=1,nS
eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) &
+ 0.5d0 * eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
+ 0.5d0 * eh_sing_rho(a,j,n)*eh_sing_rho(i,b,n)/eh_sing_Om(n) &
+ 1.5d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) &
+ 1.5d0 * eh_trip_rho(a,j,n)*eh_trip_rho(i,b,n)/eh_trip_Om(n)
end do
! do n=1,nS
! eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) &
! + 0.5d0 * eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
! + 0.5d0 * eh_sing_rho(a,j,n)*eh_sing_rho(i,b,n)/eh_sing_Om(n) &
! + 1.5d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) &
! + 1.5d0 * eh_trip_rho(a,j,n)*eh_trip_rho(i,b,n)/eh_trip_Om(n)
! end do
do n=1,nVVs
eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) &
+ ee_sing_rho(a,b,n)*ee_sing_rho(i,j,n)/ee_sing_Om(n)
end do
! do n=1,nVVs
! eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) &
! + ee_sing_rho(a,b,n)*ee_sing_rho(i,j,n)/ee_sing_Om(n)
! end do
do n=1,nOOs
eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) &
- hh_sing_rho(a,b,n)*hh_sing_rho(i,j,n)/hh_sing_Om(n)
end do
! do n=1,nOOs
! eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) &
! - hh_sing_rho(a,b,n)*hh_sing_rho(i,j,n)/hh_sing_Om(n)
! end do
do n=1,nVVt
eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) &
+ 3d0 * ee_trip_rho(a,b,n)*ee_trip_rho(i,j,n)/ee_trip_Om(n)
end do
! do n=1,nVVt
! eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) &
! + 3d0 * ee_trip_rho(a,b,n)*ee_trip_rho(i,j,n)/ee_trip_Om(n)
! end do
do n=1,nOOt
eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) &
- 3d0 * hh_trip_rho(a,b,n)*hh_trip_rho(i,j,n)/hh_trip_Om(n)
end do
! do n=1,nOOt
! eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) &
! - 3d0 * hh_trip_rho(a,b,n)*hh_trip_rho(i,j,n)/hh_trip_Om(n)
! end do
enddo
enddo

213
src/Parquet/R_eh_triplet_Gam.f90
View File

@ -1,101 +1,15 @@
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)
subroutine R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing_Phi,pp_trip_Phi,eh_trip_Gam_A)
! 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) &
+ 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) &
- 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) &
+ 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) &
+ 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) &
- 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, &
hh_sing_Om,hh_sing_rho,hh_trip_Om,hh_trip_rho, eh_trip_Gam_A)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,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)
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_trip_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: i,a,j,b
@ -118,33 +32,33 @@ subroutine R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
do b=nO+1,norb-nR
jb = jb + 1
do n=1,nS
eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) &
+ 0.5d0 * eh_sing_rho(b,a,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) &
+ 0.5d0 * eh_sing_rho(a,b,n)*eh_sing_rho(i,j,n)/eh_sing_Om(n) &
- 0.5d0 * eh_trip_rho(b,a,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n) &
- 0.5d0 * eh_trip_rho(a,b,n)*eh_trip_rho(i,j,n)/eh_trip_Om(n)
end do
! do n=1,nS
! eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) &
! + 0.5d0 * eh_sing_rho(b,a,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) &
! + 0.5d0 * eh_sing_rho(a,b,n)*eh_sing_rho(i,j,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_trip_rho(b,a,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n) &
! - 0.5d0 * eh_trip_rho(a,b,n)*eh_trip_rho(i,j,n)/eh_trip_Om(n)
! end do
do n=1,nVVs
eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) &
- ee_sing_rho(a,j,n)*ee_sing_rho(i,b,n)/ee_sing_Om(n)
end do
! do n=1,nVVs
! eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) &
! - ee_sing_rho(a,j,n)*ee_sing_rho(i,b,n)/ee_sing_Om(n)
! end do
do n=1,nOOs
eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) &
+ hh_sing_rho(a,j,n)*hh_sing_rho(i,b,n)/hh_sing_Om(n)
end do
! do n=1,nOOs
! eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) &
! + hh_sing_rho(a,j,n)*hh_sing_rho(i,b,n)/hh_sing_Om(n)
! end do
do n=1,nVVt
eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) &
+ ee_trip_rho(a,j,n)*ee_trip_rho(i,b,n)/ee_trip_Om(n)
end do
! do n=1,nVVt
! eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) &
! + ee_trip_rho(a,j,n)*ee_trip_rho(i,b,n)/ee_trip_Om(n)
! end do
do n=1,nOOt
eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) &
- hh_trip_rho(a,j,n)*hh_trip_rho(i,b,n)/hh_trip_Om(n)
end do
! do n=1,nOOt
! eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) &
! - hh_trip_rho(a,j,n)*hh_trip_rho(i,b,n)/hh_trip_Om(n)
! end do
enddo
enddo
@ -153,29 +67,18 @@ subroutine R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
end subroutine R_eh_triplet_Gamma_A
subroutine R_eh_triplet_Gamma_B(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_B)
subroutine R_eh_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,eh_sing_Phi,eh_trip_Phi,pp_sing_Phi,pp_trip_Phi,eh_trip_Gam_B)
! 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)
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_trip_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: i,a,j,b
@ -198,33 +101,33 @@ subroutine R_eh_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
do b=nO+1,norb-nR
jb = jb + 1
do n=1,nS
eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) &
+ 0.5d0 * eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
+ 0.5d0 * eh_sing_rho(a,j,n)*eh_sing_rho(i,b,n)/eh_sing_Om(n) &
- 0.5d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) &
- 0.5d0 * eh_trip_rho(a,j,n)*eh_trip_rho(i,b,n)/eh_trip_Om(n)
end do
! do n=1,nS
! eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) &
! + 0.5d0 * eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
! + 0.5d0 * eh_sing_rho(a,j,n)*eh_sing_rho(i,b,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) &
! - 0.5d0 * eh_trip_rho(a,j,n)*eh_trip_rho(i,b,n)/eh_trip_Om(n)
! end do
do n=1,nVVs
eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) &
- ee_sing_rho(a,b,n)*ee_sing_rho(i,j,n)/ee_sing_Om(n)
end do
! do n=1,nVVs
! eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) &
! - ee_sing_rho(a,b,n)*ee_sing_rho(i,j,n)/ee_sing_Om(n)
! end do
do n=1,nOOs
eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) &
+ hh_sing_rho(a,b,n)*hh_sing_rho(i,j,n)/hh_sing_Om(n)
end do
! do n=1,nOOs
! eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) &
! + hh_sing_rho(a,b,n)*hh_sing_rho(i,j,n)/hh_sing_Om(n)
! end do
do n=1,nVVt
eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) &
+ ee_trip_rho(a,b,n)*ee_trip_rho(i,j,n)/ee_trip_Om(n)
end do
! do n=1,nVVt
! eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) &
! + ee_trip_rho(a,b,n)*ee_trip_rho(i,j,n)/ee_trip_Om(n)
! end do
do n=1,nOOt
eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) &
- hh_trip_rho(a,b,n)*hh_trip_rho(i,j,n)/hh_trip_Om(n)
end do
! do n=1,nOOt
! eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) &
! - hh_trip_rho(a,b,n)*hh_trip_rho(i,j,n)/hh_trip_Om(n)
! end do
enddo
enddo

161
src/Parquet/R_pp_singlet_Gam.f90
View File

@ -1,76 +1,17 @@
subroutine R_pp_singlet_Gamma(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho,pp_sing_Gam)
subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nS,nOOs,eh_sing_Phi,eh_trip_Phi,pp_sing_Gam_D)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nR,nS
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)
! Local variables
integer :: p,q,r,s
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_sing_Gam(nOrb,nOrb,nOrb,nOrb)
! Initialization
pp_sing_Gam(:,:,:,:) = 0d0
! !$OMP PARALLEL DEFAULT(NONE) &
! !$OMP PRIVATE(i, j, ij, k, l, kl, n) &
! !$OMP SHARED(nC, nOrb, nO, nS, pp_sing_Gam_D, 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_sing_Gam(p,q,r,s) = pp_sing_Gam(p,q,r,s) &
- 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)))
end do
end do
end do
end do
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine R_pp_singlet_Gamma
subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOs,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, pp_sing_Gam_D)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOs
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)
integer,intent(in) :: nOrb,nC,nO,nS,nOOs
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: i,j,k,l
integer :: ij,kl
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_sing_Gam_D(nOOs,nOOs)
@ -93,19 +34,19 @@ subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOs,eh_sing_Om,eh_sing_rho,
do l=k,nO
kl = kl +1
do n=1,nS
pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl) &
- 0.5d0 * eh_sing_rho(k,i,n)*eh_sing_rho(j,l,n)/eh_sing_Om(n) &
- 0.5d0 * eh_sing_rho(i,k,n)*eh_sing_rho(l,j,n)/eh_sing_Om(n) &
+ 1.5d0 * eh_trip_rho(k,i,n)*eh_trip_rho(j,l,n)/eh_trip_Om(n) &
+ 1.5d0 * eh_trip_rho(i,k,n)*eh_trip_rho(l,j,n)/eh_trip_Om(n) &
- 0.5d0 * eh_sing_rho(l,i,n)*eh_sing_rho(j,k,n)/eh_sing_Om(n) &
- 0.5d0 * eh_sing_rho(i,l,n)*eh_sing_rho(k,j,n)/eh_sing_Om(n) &
+ 1.5d0 * eh_trip_rho(l,i,n)*eh_trip_rho(j,k,n)/eh_trip_Om(n) &
+ 1.5d0 * eh_trip_rho(i,l,n)*eh_trip_rho(k,j,n)/eh_trip_Om(n)
end do
! do n=1,nS
! pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl) &
! - 0.5d0 * eh_sing_rho(k,i,n)*eh_sing_rho(j,l,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_sing_rho(i,k,n)*eh_sing_rho(l,j,n)/eh_sing_Om(n) &
! + 1.5d0 * eh_trip_rho(k,i,n)*eh_trip_rho(j,l,n)/eh_trip_Om(n) &
! + 1.5d0 * eh_trip_rho(i,k,n)*eh_trip_rho(l,j,n)/eh_trip_Om(n) &
! - 0.5d0 * eh_sing_rho(l,i,n)*eh_sing_rho(j,k,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_sing_rho(i,l,n)*eh_sing_rho(k,j,n)/eh_sing_Om(n) &
! + 1.5d0 * eh_trip_rho(l,i,n)*eh_trip_rho(j,k,n)/eh_trip_Om(n) &
! + 1.5d0 * eh_trip_rho(i,l,n)*eh_trip_rho(k,j,n)/eh_trip_Om(n)
! end do
pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
! pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
end do
end do
@ -116,23 +57,20 @@ subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOs,eh_sing_Om,eh_sing_rho,
end subroutine
subroutine R_pp_singlet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVs,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, pp_sing_Gam_C)
subroutine R_pp_singlet_Gamma_C(nOrb,nO,nR,nS,nVVs,eh_sing_Phi,eh_trip_Phi,pp_sing_Gam_C)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nVVs
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)
integer,intent(in) :: nOrb,nO,nR,nS,nVVs
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: a,b,c,d
integer :: ab,cd
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_sing_Gam_C(nVVs,nVVs)
@ -155,19 +93,19 @@ subroutine R_pp_singlet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVs,eh_sing_Om,eh_sing_rho,
do d=c,nOrb - nR
cd = cd +1
do n=1,nS
pp_sing_Gam_C(ab,cd) = pp_sing_Gam_C(ab,cd) &
- 0.5d0 * eh_sing_rho(c,a,n)*eh_sing_rho(b,d,n)/eh_sing_Om(n) &
- 0.5d0 * eh_sing_rho(a,c,n)*eh_sing_rho(d,b,n)/eh_sing_Om(n) &
+ 1.5d0 * eh_trip_rho(c,a,n)*eh_trip_rho(b,d,n)/eh_trip_Om(n) &
+ 1.5d0 * eh_trip_rho(a,c,n)*eh_trip_rho(d,b,n)/eh_trip_Om(n) &
- 0.5d0 * eh_sing_rho(d,a,n)*eh_sing_rho(b,c,n)/eh_sing_Om(n) &
- 0.5d0 * eh_sing_rho(a,d,n)*eh_sing_rho(c,b,n)/eh_sing_Om(n) &
+ 1.5d0 * eh_trip_rho(d,a,n)*eh_trip_rho(b,c,n)/eh_trip_Om(n) &
+ 1.5d0 * eh_trip_rho(a,d,n)*eh_trip_rho(c,b,n)/eh_trip_Om(n)
end do
! do n=1,nS
! pp_sing_Gam_C(ab,cd) = pp_sing_Gam_C(ab,cd) &
! - 0.5d0 * eh_sing_rho(c,a,n)*eh_sing_rho(b,d,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_sing_rho(a,c,n)*eh_sing_rho(d,b,n)/eh_sing_Om(n) &
! + 1.5d0 * eh_trip_rho(c,a,n)*eh_trip_rho(b,d,n)/eh_trip_Om(n) &
! + 1.5d0 * eh_trip_rho(a,c,n)*eh_trip_rho(d,b,n)/eh_trip_Om(n) &
! - 0.5d0 * eh_sing_rho(d,a,n)*eh_sing_rho(b,c,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_sing_rho(a,d,n)*eh_sing_rho(c,b,n)/eh_sing_Om(n) &
! + 1.5d0 * eh_trip_rho(d,a,n)*eh_trip_rho(b,c,n)/eh_trip_Om(n) &
! + 1.5d0 * eh_trip_rho(a,d,n)*eh_trip_rho(c,b,n)/eh_trip_Om(n)
! end do
pp_sing_Gam_C(ab,cd) = pp_sing_Gam_C(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
! pp_sing_Gam_C(ab,cd) = pp_sing_Gam_C(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
end do
end do
@ -178,23 +116,20 @@ subroutine R_pp_singlet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVs,eh_sing_Om,eh_sing_rho,
end subroutine
subroutine R_pp_singlet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho,pp_sing_Gam_B)
subroutine R_pp_singlet_Gamma_B(nOrb,nC,nO,nR,nS,nOOs,nVVs,eh_sing_Phi,eh_trip_Phi,pp_sing_Gam_B)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOs,nVVs
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)
integer,intent(in) :: nOrb,nC,nO,nR,nS,nOOs,nVVs
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: a,b,i,j
integer :: ab,ij
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_sing_Gam_B(nVVs,nOOs)
@ -217,19 +152,19 @@ subroutine R_pp_singlet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,eh_sing_Om,eh_sing
do j=i,nO
ij = ij +1
do n=1,nS
pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij) &
- 0.5d0 * eh_sing_rho(i,a,n)*eh_sing_rho(b,j,n)/eh_sing_Om(n) &
- 0.5d0 * eh_sing_rho(a,i,n)*eh_sing_rho(j,b,n)/eh_sing_Om(n) &
+ 1.5d0 * eh_trip_rho(i,a,n)*eh_trip_rho(b,j,n)/eh_trip_Om(n) &
+ 1.5d0 * eh_trip_rho(a,i,n)*eh_trip_rho(j,b,n)/eh_trip_Om(n) &
- 0.5d0 * eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
- 0.5d0 * eh_sing_rho(a,j,n)*eh_sing_rho(i,b,n)/eh_sing_Om(n) &
+ 1.5d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) &
+ 1.5d0 * eh_trip_rho(a,j,n)*eh_trip_rho(i,b,n)/eh_trip_Om(n)
end do
! do n=1,nS
! pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij) &
! - 0.5d0 * eh_sing_rho(i,a,n)*eh_sing_rho(b,j,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_sing_rho(a,i,n)*eh_sing_rho(j,b,n)/eh_sing_Om(n) &
! + 1.5d0 * eh_trip_rho(i,a,n)*eh_trip_rho(b,j,n)/eh_trip_Om(n) &
! + 1.5d0 * eh_trip_rho(a,i,n)*eh_trip_rho(j,b,n)/eh_trip_Om(n) &
! - 0.5d0 * eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_sing_rho(a,j,n)*eh_sing_rho(i,b,n)/eh_sing_Om(n) &
! + 1.5d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) &
! + 1.5d0 * eh_trip_rho(a,j,n)*eh_trip_rho(i,b,n)/eh_trip_Om(n)
! end do
pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
! pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
end do
end do

145
src/Parquet/R_pp_triplet_Gam.f90
View File

@ -1,75 +1,17 @@
subroutine R_pp_triplet_Gamma(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, pp_trip_Gam)
subroutine R_pp_triplet_Gamma_D(nOrb,nC,nO,nS,nOOt,eh_sing_Phi,eh_trip_Phi,pp_trip_Gam_D)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nR,nS
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)
! Local variables
integer :: p,q,r,s
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_trip_Gam(nOrb,nOrb,nOrb,nOrb)
! Initialization
pp_trip_Gam(:,:,:,:) = 0d0
! !$OMP PARALLEL DEFAULT(NONE) &
! !$OMP PRIVATE(i, j, ij, k, l, kl, n) &
! !$OMP SHARED(nC, nOrb, nO, nS, pp_trip_Gam_D, 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_trip_Gam(p,q,r,s) = pp_trip_Gam(p,q,r,s) &
- 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
end do
end do
end do
! !$OMP END DO
! !$OMP END PARALLEL
end subroutine
subroutine R_pp_triplet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOt,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, pp_trip_Gam_D)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOt
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)
integer,intent(in) :: nOrb,nC,nO,nS,nOOt
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: i,j,k,l
integer :: ij,kl
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_trip_Gam_D(nOOt,nOOt)
@ -94,15 +36,16 @@ 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) &
- 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)
! pp_trip_Gam_D(ij,kl) = pp_trip_Gam_D(ij,kl) &
! - 0.5d0 * eh_sing_rho(k,i,n)*eh_sing_rho(j,l,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_sing_rho(i,k,n)*eh_sing_rho(l,j,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_trip_rho(k,i,n)*eh_trip_rho(j,l,n)/eh_trip_Om(n) &
! - 0.5d0 * eh_trip_rho(i,k,n)*eh_trip_rho(l,j,n)/eh_trip_Om(n) &
! + 0.5d0 * eh_sing_rho(l,i,n)*eh_sing_rho(j,k,n)/eh_sing_Om(n) &
! + 0.5d0 * eh_sing_rho(i,l,n)*eh_sing_rho(k,j,n)/eh_sing_Om(n) &
! + 0.5d0 * eh_trip_rho(l,i,n)*eh_trip_rho(j,k,n)/eh_trip_Om(n) &
! + 0.5d0 * eh_trip_rho(i,l,n)*eh_trip_rho(k,j,n)/eh_trip_Om(n)
end do
end do
@ -114,23 +57,20 @@ subroutine R_pp_triplet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOt,eh_sing_Om,eh_sing_rho,
end subroutine
subroutine R_pp_triplet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVt,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, pp_trip_Gam_C)
subroutine R_pp_triplet_Gamma_C(nOrb,nO,nR,nS,nVVt,eh_sing_Phi,eh_trip_Phi,pp_trip_Gam_C)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,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)
integer,intent(in) :: nOrb,nO,nR,nS,nVVt
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: a,b,c,d
integer :: ab,cd
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_trip_Gam_C(nVVt,nVVt)
@ -155,15 +95,16 @@ 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) &
- 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)
! pp_trip_Gam_C(ab,cd) = pp_trip_Gam_C(ab,cd) &
! - 0.5d0 * eh_sing_rho(c,a,n)*eh_sing_rho(b,d,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_sing_rho(a,c,n)*eh_sing_rho(d,b,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_trip_rho(c,a,n)*eh_trip_rho(b,d,n)/eh_trip_Om(n) &
! - 0.5d0 * eh_trip_rho(a,c,n)*eh_trip_rho(d,b,n)/eh_trip_Om(n) &
! + 0.5d0 * eh_sing_rho(d,a,n)*eh_sing_rho(b,c,n)/eh_sing_Om(n) &
! + 0.5d0 * eh_sing_rho(a,d,n)*eh_sing_rho(c,b,n)/eh_sing_Om(n) &
! + 0.5d0 * eh_trip_rho(d,a,n)*eh_trip_rho(b,c,n)/eh_trip_Om(n) &
! + 0.5d0 * eh_trip_rho(a,d,n)*eh_trip_rho(c,b,n)/eh_trip_Om(n)
end do
end do
@ -175,23 +116,20 @@ subroutine R_pp_triplet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVt,eh_sing_Om,eh_sing_rho,
end subroutine
subroutine R_pp_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOt,nVVt,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho,pp_trip_Gam_B)
subroutine R_pp_triplet_Gamma_B(nOrb,nC,nO,nR,nS,nOOt,nVVt,eh_sing_Phi,eh_trip_Phi,pp_trip_Gam_B)
! Compute irreducible vertex in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,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)
integer,intent(in) :: nOrb,nC,nO,nR,nS,nOOt,nVVt
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: a,b,i,j
integer :: ab,ij
integer :: n
double precision,external :: Kronecker_delta
! Output variables
double precision, intent(out) :: pp_trip_Gam_B(nVVt,nOOt)
@ -216,15 +154,16 @@ 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) &
- 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)
! pp_trip_Gam_B(ab,ij) = pp_trip_Gam_B(ab,ij) &
! - 0.5d0 * eh_sing_rho(i,a,n)*eh_sing_rho(b,j,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_sing_rho(a,i,n)*eh_sing_rho(j,b,n)/eh_sing_Om(n) &
! - 0.5d0 * eh_trip_rho(i,a,n)*eh_trip_rho(b,j,n)/eh_trip_Om(n) &
! - 0.5d0 * eh_trip_rho(a,i,n)*eh_trip_rho(j,b,n)/eh_trip_Om(n) &
! + 0.5d0 * eh_sing_rho(j,a,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) &
! + 0.5d0 * eh_sing_rho(a,j,n)*eh_sing_rho(i,b,n)/eh_sing_Om(n) &
! + 0.5d0 * eh_trip_rho(j,a,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) &
! + 0.5d0 * eh_trip_rho(a,j,n)*eh_trip_rho(i,b,n)/eh_trip_Om(n)
end do
end do

254
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,XmY,rho)
subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_trip_Phi,pp_sing_Phi,pp_trip_Phi,XpY,XmY,rho)
! Compute excitation densities
implicit none
@ -6,7 +6,10 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,XpY,X
! 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) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_trip_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: XpY(nS,nS),XmY(nS,nS)
! Local variables
@ -24,24 +27,28 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,XpY,X
! !$OMP DO
do q=nC+1,nOrb-nR
do p=nC+1,nOrb-nR
jb = 0
do j=nC+1,nO
do b=nO+1,nOrb-nR
jb = jb + 1
do ia=1,nS
X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb))
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
! X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb))
! Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
rho(p,q,ia) = rho(p,q,ia) &
+ (2d0*ERI(p,j,q,b) - ERI(p,j,b,q))*X &
+ (2d0*ERI(p,b,q,j) - ERI(p,b,j,q))*Y &
+ 1d0*eh_sing_Gam(p,j,q,b)*X &
+ 1d0*eh_sing_Gam(p,b,q,j)*Y
! rho(p,q,ia) = rho(p,q,ia) &
! + (2d0*ERI(p,j,q,b) - ERI(p,j,b,q))*X &
! + (2d0*ERI(p,b,q,j) - ERI(p,b,j,q))*Y &
! + 1d0*eh_sing_Gam(p,j,q,b)*X &
! + 1d0*eh_sing_Gam(p,b,q,j)*Y
end do
end do
end do
end do
end do
! !$OMP END DO
@ -49,7 +56,7 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,XpY,X
end subroutine R_eh_singlet_screened_integral
subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,XpY,XmY,rho)
subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_trip_Phi,pp_sing_Phi,pp_trip_Phi,XpY,XmY,rho)
! Compute excitation densities
implicit none
@ -57,7 +64,10 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,XpY,X
! 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) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_trip_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: XpY(nS,nS),XmY(nS,nS)
! Local variables
@ -81,14 +91,15 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,XpY,X
jb = jb + 1
do ia=1,nS
X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb))
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
! X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb))
! Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
! rho(p,q,ia) = rho(p,q,ia) &
! - ERI(p,j,b,q)*X &
! - ERI(p,b,j,q)*Y &
! + 1d0*eh_trip_Gam(p,j,q,b)*X &
! + 1d0*eh_trip_Gam(p,b,q,j)*Y
rho(p,q,ia) = rho(p,q,ia) &
- ERI(p,j,b,q)*X &
- ERI(p,b,j,q)*Y &
+ 1d0*eh_trip_Gam(p,j,q,b)*X &
+ 1d0*eh_trip_Gam(p,b,q,j)*Y
end do
end do
end do
@ -100,7 +111,7 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_trip_Gam,XpY,X
end subroutine R_eh_triplet_screened_integral
subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_sing_Gam,X1,Y1,rho1,X2,Y2,rho2)
subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,eh_trip_Phi,X1,Y1,rho1,X2,Y2,rho2)
! Compute excitation densities in the singlet pp channel
@ -109,11 +120,11 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_sing_G
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR
integer,intent(in) :: nOrb,nC,nO,nR
integer,intent(in) :: nOO,nVV
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_sing_Gam(nOrb,nOrb,nOrb,nOrb)
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: X1(nVV,nVV)
double precision,intent(in) :: Y1(nOO,nVV)
double precision,intent(in) :: X2(nVV,nOO)
@ -146,59 +157,62 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_sing_G
do q=nC+1,nOrb-nR
do p=nC+1,nOrb-nR
ab = 0
do a=nO+1,nOrb-nR
do b=a,nOrb-nR
ab = ab + 1
! ab = 0
! do a=nO+1,nOrb-nR
! do b=a,nOrb-nR
! ab = ab + 1
cd = 0
do c=nO+1,nOrb-nR
do d=c,nOrb-nR
cd = cd + 1
rho1(p,q,ab) = rho1(p,q,ab) &
+ (ERI(p,q,c,d) + ERI(p,q,d,c) + 1d0*pp_sing_Gam(p,q,c,d))*X1(cd,ab) &
/sqrt(1d0 + Kronecker_delta(c,d))
end do
end do
! cd = 0
! do c=nO+1,nOrb-nR
! do d=c,nOrb-nR
! cd = cd + 1
! rho1(p,q,ab) = rho1(p,q,ab) &
! + (ERI(p,q,c,d) + ERI(p,q,d,c) + 1d0*pp_sing_Gam(p,q,c,d))*X1(cd,ab) &
! /sqrt(1d0 + Kronecker_delta(c,d))
! end do
! end do
kl = 0
do k=nC+1,nO
do l=k,nO
kl = kl + 1
rho1(p,q,ab) = rho1(p,q,ab) &
+ (ERI(p,q,k,l) + ERI(p,q,l,k) + 1d0*pp_sing_Gam(p,q,k,l))*Y1(kl,ab) &
/sqrt(1d0 + Kronecker_delta(k,l))
end do
end do
! kl = 0
! do k=nC+1,nO
! do l=k,nO
! kl = kl + 1
! rho1(p,q,ab) = rho1(p,q,ab) &
! + (ERI(p,q,k,l) + ERI(p,q,l,k) + 1d0*pp_sing_Gam(p,q,k,l))*Y1(kl,ab) &
! /sqrt(1d0 + Kronecker_delta(k,l))
! end do
! end do
end do
end do
! end do
! end do
ij = 0
do i=nC+1,nO
do j=i,nO
ij = ij + 1
cd = 0
do c=nO+1,nOrb-nR
do d=c,nOrb-nR
cd = cd + 1
rho2(p,q,ij) = rho2(p,q,ij) &
+ (ERI(p,q,c,d) + ERI(p,q,d,c) + 1d0*pp_sing_Gam(p,q,c,d))*X2(cd,ij) &
/sqrt(1d0 + Kronecker_delta(c,d))
end do
end do
! ij = 0
! do i=nC+1,nO
! do j=i,nO
! ij = ij + 1
kl = 0
do k=nC+1,nO
do l=k,nO
kl = kl + 1
rho2(p,q,ij) = rho2(p,q,ij) &
+ (ERI(p,q,k,l) + ERI(p,q,l,k) + 1d0*pp_sing_Gam(p,q,k,l))*Y2(kl,ij) &
/sqrt(1d0 + Kronecker_delta(k,l))
end do
end do
end do
end do
! cd = 0
! do c=nO+1,nOrb-nR
! do d=c,nOrb-nR
! cd = cd + 1
! rho2(p,q,ij) = rho2(p,q,ij) &
! + (ERI(p,q,c,d) + ERI(p,q,d,c) + 1d0*pp_sing_Gam(p,q,c,d))*X2(cd,ij) &
! /sqrt(1d0 + Kronecker_delta(c,d))
! end do
! end do
! kl = 0
! do k=nC+1,nO
! do l=k,nO
! kl = kl + 1
! rho2(p,q,ij) = rho2(p,q,ij) &
! + (ERI(p,q,k,l) + ERI(p,q,l,k) + 1d0*pp_sing_Gam(p,q,k,l))*Y2(kl,ij) &
! /sqrt(1d0 + Kronecker_delta(k,l))
! end do
! end do
! end do
! end do
end do
end do
! !$OMP END DO
@ -209,17 +223,17 @@ end subroutine R_pp_singlet_screened_integral
subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_trip_Gam,X1,Y1,rho1,X2,Y2,rho2)
subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,eh_trip_Phi,X1,Y1,rho1,X2,Y2,rho2)
! Compute excitation densities in the triplet pp channel
implicit none
! Input variables
integer,intent(in) :: nOrb,nC,nO,nV,nR
integer,intent(in) :: nOO
integer,intent(in) :: nVV
integer,intent(in) :: nOrb,nC,nO,nR
integer,intent(in) :: nOO,nVV
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: pp_trip_Gam(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_sing_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: X1(nVV,nVV)
double precision,intent(in) :: Y1(nOO,nVV)
double precision,intent(in) :: X2(nVV,nOO)
@ -252,59 +266,61 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_trip_G
do p = nC+1, nOrb-nR
ab = 0
do a = nO+1, nOrb-nR
do b = a+1, nOrb-nR
ab = ab + 1
! do a = nO+1, nOrb-nR
! do b = a+1, nOrb-nR
! ab = ab + 1
cd = 0
do c = nO+1, nOrb-nR
do d = c+1, nOrb-nR
cd = cd + 1
! cd = 0
! do c = nO+1, nOrb-nR
! do d = c+1, nOrb-nR
! cd = cd + 1
rho1(p,q,ab) = rho1(p,q,ab) &
+ (ERI(p,q,c,d) - ERI(p,q,d,c) + 1d0*pp_trip_Gam(p,q,c,d))*X1(cd,ab)
end do ! d
end do ! c
! rho1(p,q,ab) = rho1(p,q,ab) &
! + (ERI(p,q,c,d) - ERI(p,q,d,c) + 1d0*pp_trip_Gam(p,q,c,d))*X1(cd,ab)
! end do ! d
! end do ! c
kl = 0
do k = nC+1, nO
do l = k+1, nO
! kl = 0
! do k = nC+1, nO
! do l = k+1, nO
kl = kl + 1
! kl = kl + 1
rho1(p,q,ab) = rho1(p,q,ab) &
+ (ERI(p,q,k,l) - ERI(p,q,l,k) + 1d0*pp_trip_Gam(p,q,k,l))*Y1(kl,ab)
end do ! l
end do ! k
end do ! b
end do ! a
! rho1(p,q,ab) = rho1(p,q,ab) &
! + (ERI(p,q,k,l) - ERI(p,q,l,k) + 1d0*pp_trip_Gam(p,q,k,l))*Y1(kl,ab)
! end do ! l
! end do ! k
! end do ! b
! end do ! a
ij = 0
do i = nC+1, nO
do j = i+1, nO
ij = ij + 1
! ij = 0
! do i = nC+1, nO
! do j = i+1, nO
! ij = ij + 1
cd = 0
do c = nO+1, nOrb-nR
do d = c+1, nOrb-nR
cd = cd + 1
! cd = 0
! do c = nO+1, nOrb-nR
! do d = c+1, nOrb-nR
! cd = cd + 1
rho2(p,q,ij) = rho2(p,q,ij) &
+ (ERI(p,q,c,d) - ERI(p,q,d,c) + 1d0*pp_trip_Gam(p,q,c,d))*X2(cd,ij)
end do ! d
end do ! c
! rho2(p,q,ij) = rho2(p,q,ij) &
! + (ERI(p,q,c,d) - ERI(p,q,d,c) + 1d0*pp_trip_Gam(p,q,c,d))*X2(cd,ij)
! end do ! d
! end do ! c
kl = 0
do k = nC+1, nO
do l = k+1, nO
kl = kl + 1
! kl = 0
! do k = nC+1, nO
! do l = k+1, nO
! kl = kl + 1
rho2(p,q,ij) = rho2(p,q,ij) &
+ (ERI(p,q,k,l) - ERI(p,q,l,k) + 1d0*pp_trip_Gam(p,q,k,l))*Y2(kl,ij)
end do ! l
end do ! k
end do ! j
end do ! i
! rho2(p,q,ij) = rho2(p,q,ij) &
! + (ERI(p,q,k,l) - ERI(p,q,l,k) + 1d0*pp_trip_Gam(p,q,k,l))*Y2(kl,ij)
! end do ! l
! end do ! k
! end do ! j
! end do ! i
end do ! p
end do ! q
! !$OMP END DO