4
1
mirror of https://github.com/pfloos/quack synced 2024-06-18 11:15:30 +02:00

fix bugs in print in UKS

This commit is contained in:
Pierre-Francois Loos 2022-01-08 13:59:45 +01:00
parent 41d160d1d6
commit d790d95c4d
5 changed files with 47 additions and 28 deletions

View File

@ -1,7 +1,7 @@
# RHF UHF KS MOM
T F F F
# MP2* MP3 MP2-F12
T F F
F F F
# CCD pCCD DCD CCSD CCSD(T)
F F F F F
# drCCD rCCD crCCD lCCD

View File

@ -80,10 +80,10 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
! Dimensions of the pp-RPA linear reponse matrices
! nOOs = nO*(nO + 1)/2
! nVVs = nV*(nV + 1)/2
nOOs = nO*nO
nVVs = nV*nV
nOOs = nO*(nO + 1)/2
nVVs = nV*(nV + 1)/2
! nOOs = nO*nO
! nVVs = nV*nV
nOOt = nO*(nO - 1)/2
nVVt = nV*(nV - 1)/2
@ -103,8 +103,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
!----------------------------------------------
ispin = 1
! iblock = 1
iblock = 3
iblock = 1
! iblock = 3
! Compute linear response
@ -113,8 +113,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
! EcRPA(ispin) = 1d0*EcRPA(ispin)
! call print_excitation('pp-RPA (N+2)',iblock,nVVs,Omega1s(:))
! call print_excitation('pp-RPA (N-2)',iblock,nOOs,Omega2s(:))
call print_excitation('pp-RPA (N+2)',iblock,nVVs,Omega1s(:))
call print_excitation('pp-RPA (N-2)',iblock,nOOs,Omega2s(:))
!----------------------------------------------
! alpha-alpha block
@ -131,8 +131,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
! EcRPA(ispin) = 2d0*EcRPA(ispin)
! EcRPA(ispin) = 3d0*EcRPA(ispin)
! call print_excitation('pp-RPA (N+2)',iblock,nVVt,Omega1t(:))
! call print_excitation('pp-RPA (N-2)',iblock,nOOt,Omega2t(:))
call print_excitation('pp-RPA (N+2)',iblock,nVVt,Omega1t(:))
call print_excitation('pp-RPA (N-2)',iblock,nOOt,Omega2t(:))
!----------------------------------------------
! Compute T-matrix version of the self-energy
@ -142,12 +142,12 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
SigT(:) = 0d0
Z(:) = 0d0
! iblock = 1
iblock = 3
iblock = 1
! iblock = 3
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT)
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,sqrt(1d0)*Omega1s,rho1s,sqrt(1d0)*Omega2s,rho2s,EcGM,SigT)
call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,Z)
@ -155,7 +155,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT)
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,sqrt(1.5d0)*rho1t,Omega2t,sqrt(1.5d0)*rho2t,EcGM,SigT)
call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,Z)
@ -189,8 +189,9 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
! Compute the ppRPA correlation energy
ispin = 1
iblock = 3
! iblock = 1
! iblock = 3
iblock = 1
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eG0T0,ERI_MO, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
ispin = 2

View File

@ -22,8 +22,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
! Local variables
integer :: k,l
integer :: c,d
integer :: i,j,k,l
integer :: a,b,c,d
integer :: p,q
integer :: ab,cd,ij,kl
double precision,external :: Kronecker_delta
@ -47,14 +47,20 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do ab=1,nVV
! do ab=1,nVV
ab = 0
do a=nO+1,nBas-nR
do b=a,nBas-nR
ab = ab + 1
cd = 0
do c=nO+1,nBas-nR
do d=c,nBas-nR
cd = cd + 1
rho1(p,q,ab) = rho1(p,q,ab) &
+ (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(c,d)))
+ ERI(p,q,c,d)*X1(cd,ab)
! + ERI(p,q,c,d)*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
! + (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
end do
end do
@ -63,20 +69,29 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
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))*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(k,l)))
+ ERI(p,q,k,l)*Y1(kl,ab)
! + ERI(p,q,k,l)*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(k,l)))
! + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(k,l)))
end do
end do
end do
end do
do ij=1,nOO
! do ij=1,nOO
ij = 0
do i=nC+1,nO
do j=i,nO
ij = ij + 1
cd = 0
do c=nO+1,nBas-nR
do d=c,nBas-nR
cd = cd + 1
rho2(p,q,ij) = rho2(p,q,ij) &
+ (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(c,d)))
+ ERI(p,q,c,d)*X2(cd,ij)
! + ERI(p,q,c,d)*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(c,d)))
! + (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(c,d)))
end do
end do
@ -85,11 +100,14 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
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))*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(k,l)))
+ ERI(p,q,k,l)*Y2(kl,ij)
! + ERI(p,q,k,l)*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
! + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
end do
end do
end do
end do
end do
end do

View File

@ -66,7 +66,7 @@ subroutine print_UKS(nBas,nEns,occnum,Ov,wEns,eps,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipol
HOMOb = -huge(0d0)
if(iHOMOb > 0) HOMOb = eps(iHOMOb,2)
LUMOb = +huge(0d0)
if(iLUMOb <= nBas) LUMOb = eps(iLUMOb,1)
if(iLUMOb <= nBas) LUMOb = eps(iLUMOb,2)
HOMO = max(HOMOa,HOMOb)
LUMO = min(LUMOa,LUMOb)

View File

@ -149,12 +149,12 @@ subroutine print_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,LZH,LZx,LZc,
write(*,'(A40,F16.10,A3)') ' H Levy-Zahariev shift: ',sum(LZH(:)),' au'
write(*,'(A40,F16.10,A3)') ' x Levy-Zahariev shift: ',sum(LZx(:)),' au'
write(*,'(A40,F16.10,A3)') ' c Levy-Zahariev shift: ',sum(LZc(:)),' au'
write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',sum(LZH(:))+sum(LZx(:))+sum(LZx(:)),' au'
write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',sum(LZH(:))+sum(LZx(:))+sum(LZc(:)),' au'
write(*,*)
write(*,'(A40,F16.10,A3)') ' H Levy-Zahariev shift: ',sum(LZH(:))*HaToeV,' eV'
write(*,'(A40,F16.10,A3)') ' x Levy-Zahariev shift: ',sum(LZx(:))*HaToeV,' eV'
write(*,'(A40,F16.10,A3)') ' c Levy-Zahariev shift: ',sum(LZc(:))*HaToeV,' eV'
write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',(sum(LZH(:))+sum(LZx(:))+sum(LZx(:)))*HaToeV,' eV'
write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',(sum(LZH(:))+sum(LZx(:))+sum(LZc(:)))*HaToeV,' eV'
write(*,'(A60)') '-------------------------------------------------'
write(*,*)