mirror of
https://github.com/pfloos/quack
synced 2025-04-04 07:51:36 +02:00
spin-orbital implementation of self-energy
This commit is contained in:
parent
1c279e528e
commit
a45f109a96
@ -14,7 +14,7 @@ subroutine print_excitation_energies(method,manifold,nS,Om)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer,parameter :: maxS = 25
|
||||
integer,parameter :: maxS = 50
|
||||
integer :: m
|
||||
|
||||
write(*,*)
|
||||
|
||||
@ -11,7 +11,8 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
logical :: linearize = .true.
|
||||
logical :: print_phLR = .false.
|
||||
logical :: print_ppLR = .false.
|
||||
|
||||
double precision :: eta = 100d0
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: max_it_1b,max_it_2b
|
||||
@ -77,7 +78,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
|
||||
! DIIS parameters
|
||||
|
||||
max_diis = 1
|
||||
max_diis = 2
|
||||
n_diis = 0
|
||||
rcond = 1d0
|
||||
|
||||
@ -126,21 +127,18 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
|
||||
n_it_1b = 0
|
||||
err_1b = 1d0
|
||||
|
||||
n_it_2b = 0
|
||||
err_2b = 1d0
|
||||
|
||||
|
||||
eQP(:) = eHF(:)
|
||||
eOld(:) = eHF(:)
|
||||
|
||||
|
||||
eh_rho(:,:,:) = 0d0
|
||||
ee_rho(:,:,:) = 0d0
|
||||
hh_rho(:,:,:) = 0d0
|
||||
|
||||
|
||||
old_eh_Om(:) = 0d0
|
||||
old_ee_Om(:) = 0d0
|
||||
old_hh_Om(:) = 0d0
|
||||
|
||||
|
||||
old_eh_Phi(:,:,:,:) = 0d0
|
||||
old_pp_Phi(:,:,:,:) = 0d0
|
||||
|
||||
@ -159,6 +157,11 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
write(*,*)'====================================='
|
||||
write(*,*)
|
||||
|
||||
! Initialization
|
||||
|
||||
n_it_2b = 0
|
||||
err_2b = 1d0
|
||||
|
||||
!-----------------------------------------!
|
||||
! Main loop for two-body self-consistency !
|
||||
!-----------------------------------------!
|
||||
@ -186,7 +189,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
|
||||
call wall_time(start_t)
|
||||
|
||||
call phGLR_A(.false.,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
|
||||
call phGLR_A(.false.,nOrb,nC,nO,nV,nR,nS,1d0,eOld,ERI,Aph)
|
||||
if(.not.TDA) call phGLR_B(.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
|
||||
if(n_it_2b == 1) then
|
||||
@ -236,8 +239,8 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
|
||||
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)
|
||||
call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eOld,ERI,Cpp)
|
||||
call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eOld,ERI,Dpp)
|
||||
|
||||
if(n_it_2b == 1) then
|
||||
|
||||
@ -383,7 +386,6 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
call DIIS_extrapolation(rcond,2*nOrb**4,2*nOrb**4,n_diis,err_diis,Phi_diis,err,Phi)
|
||||
print*,rcond
|
||||
|
||||
end if
|
||||
|
||||
@ -439,7 +441,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
write(*,*)' Two-body convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
stop
|
||||
!stop
|
||||
|
||||
else
|
||||
|
||||
@ -459,15 +461,19 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
|
||||
write(*,*) 'Building self-energy'
|
||||
|
||||
|
||||
call wall_time(start_t)
|
||||
!call G_irred_Parquet_self_energy(nOrb,nC,nO,nV,nR,eOld,EcGM,SigC,Z)
|
||||
call G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eOld,ERI, &
|
||||
eh_rho,old_eh_Om,ee_rho,old_ee_Om,hh_rho,old_hh_Om,EcGM,SigC,Z)
|
||||
call wall_time(end_t)
|
||||
t = end_t - start_t
|
||||
write(*,'(A50,1X,F9.3,A8)') 'Wall time for self energy =',t,' seconds'
|
||||
write(*,*)
|
||||
|
||||
eQPlin(:) = eHF(:) !+ Z(:)*SigC(:)
|
||||
|
||||
eQPlin(:) = eHF(:) + Z(:)*SigC(:)
|
||||
|
||||
call print_RG0F2(nOrb,nO,eHF,SigC,eQPlin,Z,0d0,0d0,0d0)
|
||||
|
||||
! Solve the quasi-particle equation
|
||||
|
||||
if(linearize) then
|
||||
|
||||
@ -23,7 +23,7 @@ subroutine G_eh_Phi(nOrb,nC,nR,nS,eh_Om,eh_rho,eh_Phi)
|
||||
do q = nC+1, nOrb-nR
|
||||
do p = nC+1, nOrb-nR
|
||||
|
||||
do n=1,nS
|
||||
do n=1,nS
|
||||
eh_Phi(p,q,r,s) = eh_Phi(p,q,r,s) &
|
||||
- eh_rho(r,p,n)*eh_rho(q,s,n)/eh_Om(n) &
|
||||
- eh_rho(p,r,n)*eh_rho(s,q,n)/eh_Om(n)
|
||||
|
||||
297
src/Parquet/G_irred_Parquet_self_energy.f90
Normal file
297
src/Parquet/G_irred_Parquet_self_energy.f90
Normal file
@ -0,0 +1,297 @@
|
||||
subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
|
||||
eh_rho,eh_Om,ee_rho,ee_Om,hh_rho,hh_Om,EcGM,SigC,Z)
|
||||
|
||||
! Compute correlation part of the self-energy coming from irreducible vertices contribution
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC, nO, nV, nR
|
||||
integer,intent(in) :: nS, nOO, nVV
|
||||
double precision,intent(in) :: eQP(nOrb)
|
||||
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
double precision,intent(in) :: eh_rho(nOrb,nOrb,nS)
|
||||
double precision,intent(in) :: eh_Om(nS)
|
||||
double precision,intent(in) :: ee_rho(nOrb,nOrb,nVV)
|
||||
double precision,intent(in) :: ee_Om(nVV)
|
||||
double precision,intent(in) :: hh_rho(nOrb,nOrb,nOO)
|
||||
double precision,intent(in) :: hh_Om(nOO)
|
||||
|
||||
! Local variables
|
||||
integer :: i,j,k,a,b,c
|
||||
integer :: p,n
|
||||
double precision :: eps,dem1,dem2,reg,reg1,reg2
|
||||
double precision :: num
|
||||
|
||||
! Output variables
|
||||
double precision,intent(out) :: SigC(nOrb)
|
||||
double precision,intent(out) :: Z(nOrb)
|
||||
double precision,intent(out) :: EcGM
|
||||
|
||||
! Initialize
|
||||
|
||||
SigC(:) = 0d0
|
||||
Z(:) = 0d0
|
||||
EcGM = 0d0
|
||||
|
||||
!-----------------------------!
|
||||
! GF2 part of the self-energy !
|
||||
!-----------------------------!
|
||||
do p=nC+1,nOrb-nR
|
||||
! 2h1p sum
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=nO+1,nOrb-nR
|
||||
|
||||
eps = eQP(p) + eQP(a) - eQP(i) - eQP(j)
|
||||
reg = (1d0 - exp(- 2d0 * eta * eps * eps))
|
||||
num = 0.5d0*(ERI(p,a,i,j) - ERI(p,a,j,i))**2
|
||||
|
||||
SigC(p) = SigC(p) + num*reg/eps
|
||||
Z(p) = Z(p) - num*reg/eps**2
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
! 2p1h sum
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nOrb-nR
|
||||
do b=nO+1,nOrb-nR
|
||||
|
||||
eps = eQP(p) + eQP(i) - eQP(a) - eQP(b)
|
||||
reg = (1d0 - exp(- 2d0 * eta * eps * eps))
|
||||
num = 0.5d0*(ERI(p,i,a,b) - ERI(p,i,b,a))**2
|
||||
|
||||
SigC(p) = SigC(p) + num*reg/eps
|
||||
Z(p) = Z(p) - num*reg/eps**2
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!-----------------------------!
|
||||
! eh part of the self-energy !
|
||||
!-----------------------------!
|
||||
|
||||
do p=nC+1,nOrb-nR
|
||||
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nOrb-nR
|
||||
|
||||
do n=1,nS
|
||||
!3h2p
|
||||
do j=nC+1,nO
|
||||
num = ERI(p,a,j,i) * &
|
||||
(eh_rho(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,p,n))
|
||||
|
||||
dem1 = eQP(a) - eQP(i) - eh_Om(n)
|
||||
dem2 = eQP(p) - eQP(j) + eh_Om(n)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
num = ERI(p,a,j,i) * &
|
||||
(eh_rho(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,p,n))
|
||||
|
||||
dem1 = eQP(a) - eQP(i) - eh_Om(n)
|
||||
dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
num = ERI(p,i,j,a) * &
|
||||
(eh_rho(j,p,n) * eh_rho(a,i,n) - eh_rho(j,i,n) * eh_rho(a,p,n))
|
||||
|
||||
dem1 = eQP(a) - eQP(i) + eh_Om(n)
|
||||
dem2 = eQP(p) - eQP(j) + eh_Om(n)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
num = ERI(p,a,j,i) * &
|
||||
(eh_rho(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,p,n))
|
||||
|
||||
dem1 = eQP(a) - eQP(i) + eh_Om(n)
|
||||
dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
end do ! j
|
||||
!3p2h
|
||||
do b=nO+1,nOrb-nR
|
||||
num = ERI(p,a,b,i) * &
|
||||
(eh_rho(b,p,n) * eh_rho(i,a,n) - eh_rho(b,a,n) * eh_rho(i,p,n))
|
||||
|
||||
dem1 = eQP(a) - eQP(i) + eh_Om(n)
|
||||
dem2 = eQP(p) - eQP(b) - eh_Om(n)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
num = ERI(p,i,b,a) * &
|
||||
(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,p,n))
|
||||
|
||||
dem1 = eQP(a) - eQP(i) + eh_Om(n)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
num = ERI(p,i,b,a) * &
|
||||
(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,p,n))
|
||||
|
||||
dem1 = eQP(a) - eQP(i) - eh_Om(n)
|
||||
dem2 = eQP(p) - eQP(b) - eh_Om(n)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
num = ERI(p,i,b,a) * &
|
||||
(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,p,n))
|
||||
|
||||
dem1 = eQP(a) - eQP(i) - eh_Om(n)
|
||||
dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
end do ! b
|
||||
|
||||
end do ! n
|
||||
|
||||
end do ! a
|
||||
end do ! i
|
||||
|
||||
end do ! p
|
||||
|
||||
!-----------------------------!
|
||||
! pp part of the self-energy !
|
||||
!-----------------------------!
|
||||
|
||||
do p=nC+1,nOrb-nR
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do n=1,nVV
|
||||
! 4h1p
|
||||
do k=nC+1,nO
|
||||
num = ERI(p,k,i,j) * ee_rho(i,j,n) * ee_rho(p,k,n)
|
||||
dem1 = ee_Om(n) - eQP(i) - eQP(j)
|
||||
dem2 = eQP(p) + eQP(k) - ee_Om(n)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
end do ! k
|
||||
! 3h2p
|
||||
do c=nO+1,nOrb-nR
|
||||
|
||||
num = ERI(p,c,i,j) * ee_rho(i,j,n) * ee_rho(p,c,n)
|
||||
dem1 = ee_Om(n) - eQP(i) - eQP(j)
|
||||
dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
end do ! a
|
||||
end do ! n
|
||||
do n=1,nOO
|
||||
! 3h2p
|
||||
do c=nO+1,nOrb-nR
|
||||
|
||||
num = ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n)
|
||||
dem1 = hh_Om(n) - eQP(i) - eQP(j)
|
||||
dem2 = eQP(p) + eQP(c) - hh_Om(n)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
num = ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n)
|
||||
dem1 = hh_Om(n) - eQP(i) - eQP(j)
|
||||
dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
end do ! c
|
||||
end do ! n
|
||||
end do ! j
|
||||
end do ! i
|
||||
|
||||
do a=nO+1,nOrb-nR
|
||||
do b=nO+1,nOrb-nR
|
||||
do n=1,nOO
|
||||
! 4p1h
|
||||
do c=nO+1,nOrb-nR
|
||||
|
||||
num = ERI(p,c,a,b) * hh_rho(a,b,n) * hh_rho(p,c,n)
|
||||
dem1 = hh_Om(n) - eQP(a) - eQP(b)
|
||||
dem2 = eQP(p) + eQP(c) - hh_Om(n)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
end do ! c
|
||||
! 3p2h
|
||||
do k=nC+1,nO
|
||||
|
||||
num = ERI(p,k,a,b) * hh_rho(a,b,n) * hh_rho(p,k,n)
|
||||
dem1 = hh_Om(n) - eQP(a) - eQP(b)
|
||||
dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
end do ! k
|
||||
end do ! n
|
||||
do n=1,nVV
|
||||
! 3p2h
|
||||
do k=nC+1,nO
|
||||
|
||||
num = ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n)
|
||||
dem1 = ee_Om(n) - eQP(a) - eQP(b)
|
||||
dem2 = eQP(p) + eQP(k) - ee_Om(n)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) - num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
num = ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n)
|
||||
dem1 = ee_Om(n) - eQP(a) - eQP(b)
|
||||
dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b)
|
||||
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
|
||||
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
|
||||
|
||||
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
|
||||
|
||||
end do ! c
|
||||
end do ! n
|
||||
end do ! b
|
||||
end do ! a
|
||||
|
||||
end do ! p
|
||||
|
||||
!-----------------------------!
|
||||
! Renormalization factor !
|
||||
!-----------------------------!
|
||||
|
||||
Z(:) = 1d0/(1d0 - Z(:))
|
||||
|
||||
end subroutine G_Parquet_self_energy
|
||||
@ -30,7 +30,7 @@ subroutine G_pp_Phi(nOrb,nC,nR,nOO,nVV,ee_Om,ee_rho,hh_Om,hh_rho,pp_Phi)
|
||||
+ 2d0 * ee_rho(p,q,n)*ee_rho(r,s,n)/ee_Om(n)
|
||||
end do
|
||||
|
||||
do n=1,nOO
|
||||
do n=1,nOO
|
||||
pp_Phi(p,q,r,s) = pp_Phi(p,q,r,s) &
|
||||
- 2d0 * hh_rho(p,q,n)*hh_rho(r,s,n)/hh_Om(n)
|
||||
end do
|
||||
|
||||
@ -37,10 +37,12 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho
|
||||
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
|
||||
|
||||
rho(p,q,ia) = rho(p,q,ia) &
|
||||
|
||||
+ (1d0*ERI(q,j,p,b) - 1d0*ERI(q,j,b,p) &
|
||||
- 1d0*eh_Phi(q,j,b,p) + 1d0*pp_Phi(q,j,p,b)) * X &
|
||||
- 0d0*eh_Phi(q,j,b,p) + 0d0*pp_Phi(q,j,p,b)) * X &
|
||||
+ (1d0*ERI(q,b,p,j) - 1d0*ERI(q,b,j,p) &
|
||||
- 1d0*eh_Phi(q,b,j,p) + 1d0*pp_Phi(q,b,p,j)) * Y
|
||||
- 0d0*eh_Phi(q,b,j,p) + 0d0*pp_Phi(q,b,p,j)) * Y
|
||||
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -107,8 +109,10 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
|
||||
do c=nO+1,nOrb-nR
|
||||
do d=c+1,nOrb-nR
|
||||
cd = cd + 1
|
||||
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ( 1d0*ERI(p,q,c,d) - 1d0*ERI(p,q,d,c) &
|
||||
+ 1d0*eh_Phi(p,q,c,d) - 1d0*eh_Phi(p,q,d,c) ) * X1(cd,ab)
|
||||
+ 0d0*eh_Phi(p,q,c,d) - 0d0*eh_Phi(p,q,d,c) ) * X1(cd,ab)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -116,8 +120,10 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
|
||||
do k=nC+1,nO
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ( 1d0*ERI(p,q,k,l) - 1d0*ERI(p,q,l,k) &
|
||||
+ 1d0*eh_Phi(p,q,k,l) - 1d0*eh_Phi(p,q,l,k) ) * Y1(kl,ab)
|
||||
+ 0d0*eh_Phi(p,q,k,l) - 0d0*eh_Phi(p,q,l,k) ) * Y1(kl,ab)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -132,8 +138,10 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
|
||||
do c=nO+1,nOrb-nR
|
||||
do d=c+1,nOrb-nR
|
||||
cd = cd + 1
|
||||
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ( 1d0*ERI(p,q,c,d) - 1d0*ERI(p,q,d,c) &
|
||||
+ 1d0*eh_Phi(p,q,c,d) - 1d0*eh_Phi(p,q,d,c) ) * X2(cd,ij)
|
||||
+ 0d0*eh_Phi(p,q,c,d) - 0d0*eh_Phi(p,q,d,c) ) * X2(cd,ij)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -141,8 +149,10 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
|
||||
do k=nC+1,nO
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ( 1d0*ERI(p,q,k,l) - 1d0*ERI(p,q,l,k) &
|
||||
+ 1d0*eh_Phi(p,q,k,l) - 1d0*eh_Phi(p,q,l,k) ) * Y2(kl,ij)
|
||||
+ 0d0*eh_Phi(p,q,k,l) - 0d0*eh_Phi(p,q,l,k) ) * Y2(kl,ij)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -39,11 +39,11 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr
|
||||
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
|
||||
|
||||
rho(p,q,ia) = rho(p,q,ia) + (2d0*ERI(q,j,p,b) - ERI(q,j,b,p) &
|
||||
- 0.5d0*eh_sing_Phi(q,j,b,p) - 1.5d0*eh_trip_Phi(q,j,b,p) &
|
||||
+ 0.5d0*pp_sing_Phi(q,j,p,b) + 1.5d0*pp_trip_Phi(q,j,p,b)) * X &
|
||||
- 0d0*0.5d0*eh_sing_Phi(q,j,b,p) - 0d0*1.5d0*eh_trip_Phi(q,j,b,p) &
|
||||
+ 0d0*0.5d0*pp_sing_Phi(q,j,p,b) + 0d0*1.5d0*pp_trip_Phi(q,j,p,b)) * X &
|
||||
+ (2d0*ERI(q,b,p,j) - ERI(q,b,j,p) &
|
||||
- 0.5d0*eh_sing_Phi(q,b,j,p) - 1.5d0*eh_trip_Phi(q,b,j,p) &
|
||||
+ 0.5d0*pp_sing_Phi(q,b,p,j) + 1.5d0*pp_trip_Phi(q,b,p,j)) * Y
|
||||
- 0d0*0.5d0*eh_sing_Phi(q,b,j,p) - 0d0*1.5d0*eh_trip_Phi(q,b,j,p) &
|
||||
+ 0d0*0.5d0*pp_sing_Phi(q,b,p,j) + 0d0*1.5d0*pp_trip_Phi(q,b,p,j)) * Y
|
||||
|
||||
end do
|
||||
|
||||
@ -98,11 +98,11 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr
|
||||
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
|
||||
|
||||
rho(p,q,ia) = rho(p,q,ia) + (- ERI(q,j,b,p) &
|
||||
- 0.5d0*eh_sing_Phi(q,j,b,p) + 0.5d0*eh_trip_Phi(q,j,b,p) &
|
||||
- 0.5d0*pp_sing_Phi(q,j,p,b) + 0.5d0*pp_trip_Phi(q,j,p,b)) * X &
|
||||
- 0d0*0.5d0*eh_sing_Phi(q,j,b,p) + 0d0*0.5d0*eh_trip_Phi(q,j,b,p) &
|
||||
- 0d0*0.5d0*pp_sing_Phi(q,j,p,b) + 0d0*0.5d0*pp_trip_Phi(q,j,p,b)) * X &
|
||||
+ (- ERI(q,b,j,p) &
|
||||
- 0.5d0*eh_sing_Phi(q,b,j,p) + 0.5d0*eh_trip_Phi(q,b,j,p) &
|
||||
- 0.5d0*pp_sing_Phi(q,b,p,j) + 0.5d0*pp_trip_Phi(q,b,p,j)) * Y
|
||||
- 0d0*0.5d0*eh_sing_Phi(q,b,j,p) + 0d0*0.5d0*eh_trip_Phi(q,b,j,p) &
|
||||
- 0d0*0.5d0*pp_sing_Phi(q,b,p,j) + 0d0*0.5d0*pp_trip_Phi(q,b,p,j)) * Y
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -174,8 +174,8 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
|
||||
cd = cd + 1
|
||||
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) + ERI(p,q,d,c) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,c,d) - 1.5d0*eh_trip_Phi(p,q,c,d) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,d,c) - 1.5d0*eh_trip_Phi(p,q,d,c))&
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) - 0d0*1.5d0*eh_trip_Phi(p,q,c,d) &
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*1.5d0*eh_trip_Phi(p,q,d,c))&
|
||||
*X1(cd,ab)/sqrt(1d0 + Kronecker_delta(c,d))
|
||||
|
||||
end do ! d
|
||||
@ -188,8 +188,8 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
|
||||
kl = kl + 1
|
||||
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) + ERI(p,q,l,k) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,k,l) - 1.5d0*eh_trip_Phi(p,q,k,l) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,l,k) - 1.5d0*eh_trip_Phi(p,q,l,k))&
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) - 0d0*1.5d0*eh_trip_Phi(p,q,k,l) &
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*1.5d0*eh_trip_Phi(p,q,l,k))&
|
||||
*Y1(kl,ab)/sqrt(1d0 + Kronecker_delta(k,l))
|
||||
end do ! l
|
||||
end do ! k
|
||||
@ -207,8 +207,8 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
|
||||
cd = cd + 1
|
||||
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) + ERI(p,q,d,c) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,c,d) - 1.5d0*eh_trip_Phi(p,q,c,d) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,d,c) - 1.5d0*eh_trip_Phi(p,q,d,c))&
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) - 0d0*1.5d0*eh_trip_Phi(p,q,c,d) &
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*1.5d0*eh_trip_Phi(p,q,d,c))&
|
||||
*X2(cd,ij)/sqrt(1d0 + Kronecker_delta(c,d))
|
||||
end do ! d
|
||||
end do ! c
|
||||
@ -219,8 +219,8 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
|
||||
kl = kl + 1
|
||||
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) + ERI(p,q,l,k) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,k,l) - 1.5d0*eh_trip_Phi(p,q,k,l) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,l,k) - 1.5d0*eh_trip_Phi(p,q,l,k))&
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) - 0d0*1.5d0*eh_trip_Phi(p,q,k,l) &
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*1.5d0*eh_trip_Phi(p,q,l,k))&
|
||||
*Y2(kl,ij)/sqrt(1d0 + Kronecker_delta(k,l))
|
||||
end do ! l
|
||||
end do ! k
|
||||
@ -291,8 +291,8 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
|
||||
cd = cd + 1
|
||||
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) - ERI(p,q,d,c) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,c,d) + 0.5d0*eh_trip_Phi(p,q,c,d) &
|
||||
- 0.5d0*eh_sing_Phi(p,q,d,c) - 0.5d0*eh_trip_Phi(p,q,d,c) )*X1(cd,ab)
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) + 0d0*0.5d0*eh_trip_Phi(p,q,c,d) &
|
||||
- 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*0.5d0*eh_trip_Phi(p,q,d,c) )*X1(cd,ab)
|
||||
|
||||
end do ! d
|
||||
end do ! c
|
||||
@ -304,8 +304,8 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
|
||||
kl = kl + 1
|
||||
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) - ERI(p,q,l,k) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,k,l) + 0.5d0*eh_trip_Phi(p,q,k,l) &
|
||||
- 0.5d0*eh_sing_Phi(p,q,l,k) - 0.5d0*eh_trip_Phi(p,q,l,k) )*Y1(kl,ab)
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) + 0d0*0.5d0*eh_trip_Phi(p,q,k,l) &
|
||||
- 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*0.5d0*eh_trip_Phi(p,q,l,k) )*Y1(kl,ab)
|
||||
end do ! l
|
||||
end do ! k
|
||||
end do ! b
|
||||
@ -322,8 +322,8 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
|
||||
cd = cd + 1
|
||||
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) - ERI(p,q,d,c) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,c,d) + 0.5d0*eh_trip_Phi(p,q,c,d) &
|
||||
- 0.5d0*eh_sing_Phi(p,q,d,c) - 0.5d0*eh_trip_Phi(p,q,d,c) )*X2(cd,ij)
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) + 0d0*0.5d0*eh_trip_Phi(p,q,c,d) &
|
||||
- 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*0.5d0*eh_trip_Phi(p,q,d,c) )*X2(cd,ij)
|
||||
end do ! d
|
||||
end do ! c
|
||||
|
||||
@ -333,8 +333,8 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
|
||||
kl = kl + 1
|
||||
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) - ERI(p,q,l,k) &
|
||||
+ 0.5d0*eh_sing_Phi(p,q,k,l) + 0.5d0*eh_trip_Phi(p,q,k,l) &
|
||||
- 0.5d0*eh_sing_Phi(p,q,l,k) - 0.5d0*eh_trip_Phi(p,q,l,k) )*Y2(kl,ij)
|
||||
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) + 0d0*0.5d0*eh_trip_Phi(p,q,k,l) &
|
||||
- 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*0.5d0*eh_trip_Phi(p,q,l,k) )*Y2(kl,ij)
|
||||
end do ! l
|
||||
end do ! k
|
||||
|
||||
|
||||
@ -379,7 +379,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,
|
||||
call wall_time(start_Parquet)
|
||||
call RParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body, &
|
||||
nOrb,nC,nO,nV,nR,nS, &
|
||||
eHF,ERI_MO)
|
||||
eGW,ERI_MO)
|
||||
call wall_time(end_Parquet)
|
||||
|
||||
t_Parquet = end_Parquet - start_Parquet
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user