mirror of
https://github.com/pfloos/quack
synced 2024-06-16 18:25:29 +02:00
fixing up qsGW Ec
This commit is contained in:
parent
a7594c270a
commit
5f55c1af4d
|
|
@ -6,14 +6,14 @@
|
|||
# GGA = 2: B88,G96,PBE
|
||||
# MGGA = 3:
|
||||
# Hybrid = 4: HF,B3,PBE
|
||||
2 B88
|
||||
1 S51
|
||||
# correlation rung:
|
||||
# Hartree = 0: H
|
||||
# LDA = 1: PW92,VWN3,VWN5,eVWN5
|
||||
# GGA = 2: LYP,PBE
|
||||
# MGGA = 3:
|
||||
# Hybrid = 4: HF,LYP,PBE
|
||||
1 PW92
|
||||
0 H
|
||||
# quadrature grid SG-n
|
||||
1
|
||||
# Number of states in ensemble (nEns)
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
# RHF UHF KS MOM
|
||||
T F F F
|
||||
F T F F
|
||||
# MP2* MP3 MP2-F12
|
||||
F F F
|
||||
# CCD DCD CCSD CCSD(T)
|
||||
|
|
@ -13,7 +13,7 @@
|
|||
# G0F2 evGF2 G0F3 evGF3
|
||||
F F F F
|
||||
# G0W0* evGW* qsGW*
|
||||
F F T
|
||||
T F T
|
||||
# G0T0 evGT qsGT
|
||||
F F F
|
||||
# MCMP2
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z,EcRPA,EqsGW,dipole)
|
||||
subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z,EcGM,EcRPA,EqsGW,dipole)
|
||||
|
||||
! Print one-electron energies and other stuff for qsGW
|
||||
|
||||
|
|
@ -7,9 +7,18 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
|
|||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas,nO,nSCF
|
||||
double precision,intent(in) :: ENuc,EcRPA,Conv,thresh
|
||||
double precision,intent(in) :: eHF(nBas),eGW(nBas),c(nBas),P(nBas,nBas)
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nSCF
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: EcGM
|
||||
double precision,intent(in) :: EcRPA
|
||||
double precision,intent(in) :: Conv
|
||||
double precision,intent(in) :: thresh
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: eGW(nBas)
|
||||
double precision,intent(in) :: c(nBas)
|
||||
double precision,intent(in) :: P(nBas,nBas)
|
||||
double precision,intent(in) :: T(nBas,nBas),V(nBas,nBas)
|
||||
double precision,intent(in) :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas)
|
||||
double precision,intent(in) :: Z(nBas),SigC(nBas,nBas)
|
||||
|
|
@ -37,9 +46,8 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
|
|||
EV = trace_matrix(nBas,matmul(P,V))
|
||||
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
|
||||
Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
|
||||
Ec = 0d0
|
||||
! Ec = -0.50d0*trace_matrix(nBas,matmul(P,SigC))
|
||||
EqsGW = ET + EV + EJ + Ex + Ec
|
||||
EqsGW = ET + EV + EJ + Ex + EcGM
|
||||
|
||||
! Dump results
|
||||
|
||||
|
|
@ -69,7 +77,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
|
|||
write(*,*)'-------------------------------------------'
|
||||
write(*,'(2X,A30,F15.6,A3)') ' qsGW total energy:',EqsGW + ENuc,' au'
|
||||
write(*,'(2X,A30,F15.6,A3)') ' qsGW exchange energy:',Ex,' au'
|
||||
! write(*,'(2X,A30,F15.6,A3)') ' qsGW correlation energy:',Ec,' au'
|
||||
write(*,'(2X,A30,F15.6,A3)') ' qsGW correlation energy:',EcGM,' au'
|
||||
write(*,'(2X,A30,F15.6,A3)') 'RPA@qsGW correlation energy:',EcRPA,' au'
|
||||
write(*,*)'-------------------------------------------'
|
||||
write(*,*)
|
||||
|
|
@ -89,7 +97,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
|
|||
write(*,'(A32,1X,F16.10,A3)') ' Two-electron energy: ',EJ + Ex,' au'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Hartree energy: ',EJ,' au'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Exchange energy: ',Ex,' au'
|
||||
! write(*,'(A32,1X,F16.10,A3)') ' Correlation energy: ',Ec,' au'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Correlation energy: ',EcGM,' au'
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EqsGW,' au'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au'
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
|
||||
ENuc,ET,EV,EJ,Ex,Ec,EcRPA,EqsGW,SigC,Z,dipole)
|
||||
ENuc,ET,EV,EJ,Ex,Ec,EcGM,EcRPA,EqsGW,SigC,Z,dipole)
|
||||
|
||||
! Print one-electron energies and other stuff for qsUGW
|
||||
|
||||
|
|
@ -17,6 +17,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
|
|||
double precision,intent(in) :: EJ(nsp)
|
||||
double precision,intent(in) :: Ex(nspin)
|
||||
double precision,intent(in) :: Ec(nsp)
|
||||
double precision,intent(in) :: EcGM(nspin)
|
||||
double precision,intent(in) :: EcRPA
|
||||
double precision,intent(in) :: EqsGW
|
||||
double precision,intent(in) :: Conv
|
||||
|
|
@ -103,6 +104,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
|
|||
write(*,'(2X,A30,F15.6,A3)') ' qsUGW total energy:',EqsGW + ENuc,' au'
|
||||
write(*,'(2X,A30,F15.6,A3)') ' qsUGW exchange energy:',sum(Ex(:)),' au'
|
||||
! write(*,'(2X,A30,F15.6,A3)') ' qsUGW correlation energy:',sum(Ec(:)),' au'
|
||||
write(*,'(2X,A30,F15.6,A3)') ' GM@qsUGW correlation energy:',sum(EcGM(:)),' au'
|
||||
write(*,'(2X,A30,F15.6,A3)') 'RPA@qsUGW correlation energy:',EcRPA,' au'
|
||||
write(*,*)'-------------------------------------------------------------------------------&
|
||||
-------------------------------------------------'
|
||||
|
|
@ -141,7 +143,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
|
|||
write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',sum(Ex(:)),' au'
|
||||
write(*,'(A40,1X,F16.10,A3)') ' Exchange a energy: ',Ex(1),' au'
|
||||
write(*,'(A40,1X,F16.10,A3)') ' Exchange b energy: ',Ex(2),' au'
|
||||
! write(*,*)
|
||||
write(*,*)
|
||||
! write(*,'(A40,1X,F16.10,A3)') ' Correlation energy: ',sum(Ec(:)),' au'
|
||||
! write(*,'(A40,1X,F16.10,A3)') ' Correlation aa energy: ',Ec(1),' au'
|
||||
! write(*,'(A40,1X,F16.10,A3)') ' Correlation ab energy: ',Ec(2),' au'
|
||||
|
|
|
|||
|
|
@ -235,7 +235,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
|
|||
! Print results
|
||||
|
||||
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
|
||||
call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigCp,Z,EcRPA,EqsGW,dipole)
|
||||
call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigCp,Z,EcGM,EcRPA,EqsGW,dipole)
|
||||
|
||||
enddo
|
||||
!------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -77,6 +77,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
|
|||
double precision :: Ex(nspin)
|
||||
double precision :: Ec(nsp)
|
||||
double precision :: EcRPA
|
||||
double precision :: EcGM(nspin)
|
||||
double precision :: EqsGW
|
||||
double precision :: EcBSE(nspin)
|
||||
double precision :: EcAC(nspin)
|
||||
|
|
@ -230,12 +231,12 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
|
|||
|
||||
if(G0W) then
|
||||
|
||||
call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC)
|
||||
call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM)
|
||||
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z)
|
||||
|
||||
else
|
||||
|
||||
call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC)
|
||||
call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM)
|
||||
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z)
|
||||
|
||||
endif
|
||||
|
|
@ -321,7 +322,8 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
|
|||
! Coulomb energy
|
||||
|
||||
EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
|
||||
EJ(2) = 1.0d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2)))
|
||||
EJ(2) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) &
|
||||
+ 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,1)))
|
||||
EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
|
||||
|
||||
! Exchange energy
|
||||
|
|
@ -348,7 +350,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
|
|||
!------------------------------------------------------------------------
|
||||
|
||||
call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
|
||||
call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,Ec,EcRPA,EqsGW,SigCp,Z,dipole)
|
||||
call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,Ec,EcGM,EcRPA,EqsGW,SigCp,Z,dipole)
|
||||
|
||||
enddo
|
||||
!------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -16,7 +16,9 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
|
|||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b,p,x,y,jb
|
||||
integer :: i,j,a,b
|
||||
integer :: p,q,r
|
||||
integer :: jb
|
||||
double precision :: eps
|
||||
|
||||
! Output variables
|
||||
|
|
@ -26,7 +28,7 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
|
|||
|
||||
! Initialize
|
||||
|
||||
SigC = 0d0
|
||||
SigC(:,:) = 0d0
|
||||
|
||||
!-----------------------------!
|
||||
! COHSEX static approximation !
|
||||
|
|
@ -36,32 +38,32 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
|
|||
|
||||
! COHSEX: SEX of the COHSEX correlation self-energy
|
||||
|
||||
do x=nC+1,nBas-nR
|
||||
do y=nC+1,nBas-nR
|
||||
do p=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
do jb=1,nS
|
||||
SigC(x,y) = SigC(x,y) + 4d0*rho(x,i,jb)*rho(y,i,jb)/Omega(jb)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
SigC(p,q) = SigC(p,q) + 4d0*rho(p,i,jb)*rho(q,i,jb)/Omega(jb)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! COHSEX: COH part of the COHSEX correlation self-energy
|
||||
|
||||
do x=nC+1,nBas-nR
|
||||
do y=nC+1,nBas-nR
|
||||
do p=nC+1,nBas-nR
|
||||
do p=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do r=nC+1,nBas-nR
|
||||
do jb=1,nS
|
||||
SigC(x,y) = SigC(x,y) - 2d0*rho(x,p,jb)*rho(y,p,jb)/Omega(jb)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
SigC(p,q) = SigC(p,q) - 2d0*rho(p,r,jb)*rho(q,r,jb)/Omega(jb)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
EcGM = 0d0
|
||||
do i=nC+1,nO
|
||||
EcGM = EcGM + 0.5d0*SigC(i,i)
|
||||
enddo
|
||||
end do
|
||||
|
||||
else
|
||||
|
||||
|
|
@ -71,30 +73,42 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
|
|||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do x=nC+1,nBas-nR
|
||||
do y=nC+1,nBas-nR
|
||||
do p=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
do jb=1,nS
|
||||
eps = e(x) - e(i) + Omega(jb)
|
||||
SigC(x,y) = SigC(x,y) + 2d0*rho(x,i,jb)*rho(y,i,jb)*eps/(eps**2 + eta**2)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
eps = e(p) - e(i) + Omega(jb)
|
||||
SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do x=nC+1,nBas-nR
|
||||
do y=nC+1,nBas-nR
|
||||
do p=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do a=nO+1,nBas-nR
|
||||
do jb=1,nS
|
||||
eps = e(x) - e(a) - Omega(jb)
|
||||
SigC(x,y) = SigC(x,y) + 2d0*rho(x,a,jb)*rho(y,a,jb)*eps/(eps**2 + eta**2)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
eps = e(p) - e(a) - Omega(jb)
|
||||
SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
endif
|
||||
! GM correlation energy
|
||||
|
||||
EcGM = 0d0
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
do jb=1,nS
|
||||
eps = e(a) - e(i) + Omega(jb)
|
||||
EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end if
|
||||
|
||||
end subroutine self_energy_correlation
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC)
|
||||
subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy
|
||||
|
||||
|
|
@ -26,10 +26,12 @@ subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega
|
|||
! Output variables
|
||||
|
||||
double precision,intent(out) :: SigC(nBas,nBas,nspin)
|
||||
double precision :: EcGM(nspin)
|
||||
|
||||
! Initialize
|
||||
|
||||
SigC(:,:,:) = 0d0
|
||||
EcGM(:) = 0d0
|
||||
|
||||
!--------------!
|
||||
! Spin-up part !
|
||||
|
|
@ -61,6 +63,17 @@ subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega
|
|||
end do
|
||||
end do
|
||||
|
||||
! GM correlation energy
|
||||
|
||||
do i=nC(1)+1,nO(1)
|
||||
do a=nO(1)+1,nBas-nR(1)
|
||||
do jb=1,nSt
|
||||
eps = e(a,1) - e(i,1) + Omega(jb)
|
||||
EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!----------------!
|
||||
! Spin-down part !
|
||||
!----------------!
|
||||
|
|
@ -91,4 +104,15 @@ subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega
|
|||
end do
|
||||
end do
|
||||
|
||||
! GM correlation energy
|
||||
|
||||
do i=nC(2)+1,nO(2)
|
||||
do a=nO(2)+1,nBas-nR(2)
|
||||
do jb=1,nSt
|
||||
eps = e(a,2) - e(i,2) + Omega(jb)
|
||||
EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine unrestricted_self_energy_correlation
|
||||
|
|
|
|||
|
|
@ -58,6 +58,8 @@ subroutine UPW92_lda_correlation_energy(nGrid,weight,rho,Ec)
|
|||
|
||||
ra = max(0d0,rho(iG,1))
|
||||
rb = max(0d0,rho(iG,2))
|
||||
r = ra + rb
|
||||
z = (ra - rb)/r
|
||||
|
||||
! alpha-alpha contribution
|
||||
|
||||
|
|
@ -74,11 +76,9 @@ subroutine UPW92_lda_correlation_energy(nGrid,weight,rho,Ec)
|
|||
|
||||
! alpha-beta contribution
|
||||
|
||||
if(ra > threshold .or. rb > threshold) then
|
||||
if(r > threshold) then
|
||||
|
||||
r = ra + rb
|
||||
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
|
||||
z = (ra - rb)/r
|
||||
|
||||
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
|
||||
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
|
||||
|
|
|
|||
|
|
@ -173,6 +173,7 @@ subroutine UPW92_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
|
|||
|
||||
decdrb = decdrb_p + decdrb_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 &
|
||||
+ (decdrb_f - decdrb_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3
|
||||
|
||||
Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdrb*r)
|
||||
|
||||
end if
|
||||
|
|
|
|||
|
|
@ -52,6 +52,8 @@ subroutine UVWN3_lda_correlation_energy(nGrid,weight,rho,Ec)
|
|||
|
||||
ra = max(0d0,rho(iG,1))
|
||||
rb = max(0d0,rho(iG,2))
|
||||
r = ra + rb
|
||||
z = (ra - rb)/r
|
||||
|
||||
! alpha-alpha contribution
|
||||
|
||||
|
|
@ -73,11 +75,9 @@ subroutine UVWN3_lda_correlation_energy(nGrid,weight,rho,Ec)
|
|||
|
||||
! alpha-beta contribution
|
||||
|
||||
if(ra > threshold .or. rb > threshold) then
|
||||
if(r > threshold) then
|
||||
|
||||
r = ra + rb
|
||||
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
|
||||
z = (ra - rb)/r
|
||||
x = sqrt(rs)
|
||||
|
||||
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
|
||||
|
|
|
|||
|
|
@ -59,21 +59,22 @@ subroutine UVWN3_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
|
|||
|
||||
ra = max(0d0,rho(iG,1))
|
||||
rb = max(0d0,rho(iG,2))
|
||||
r = ra + rb
|
||||
z = (ra - rb)/r
|
||||
|
||||
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
|
||||
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
|
||||
|
||||
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
|
||||
|
||||
|
||||
! spin-up contribution
|
||||
|
||||
if(ra > threshold) then
|
||||
|
||||
r = ra + rb
|
||||
|
||||
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
|
||||
z = (ra - rb)/r
|
||||
x = sqrt(rs)
|
||||
|
||||
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
|
||||
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
|
||||
|
||||
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
|
||||
|
||||
x_p = x*x + b_p*x + c_p
|
||||
x_f = x*x + b_f*x + c_f
|
||||
x_a = x*x + b_a*x + c_a
|
||||
|
|
@ -132,16 +133,9 @@ subroutine UVWN3_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
|
|||
|
||||
if(rb > threshold) then
|
||||
|
||||
r = ra + rb
|
||||
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
|
||||
z = (ra - rb)/r
|
||||
x = sqrt(rs)
|
||||
|
||||
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
|
||||
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
|
||||
|
||||
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
|
||||
|
||||
|
||||
x_p = x*x + b_p*x + c_p
|
||||
x_f = x*x + b_f*x + c_f
|
||||
x_a = x*x + b_a*x + c_a
|
||||
|
|
|
|||
|
|
@ -52,6 +52,8 @@ subroutine UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec)
|
|||
|
||||
ra = max(0d0,rho(iG,1))
|
||||
rb = max(0d0,rho(iG,2))
|
||||
r = ra + rb
|
||||
z = (ra - rb)/r
|
||||
|
||||
! alpha-alpha contribution
|
||||
|
||||
|
|
@ -73,11 +75,9 @@ subroutine UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec)
|
|||
|
||||
! alpha-beta contribution
|
||||
|
||||
if(ra > threshold .or. rb > threshold) then
|
||||
if(r > threshold) then
|
||||
|
||||
r = ra + rb
|
||||
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
|
||||
z = (ra - rb)/r
|
||||
x = sqrt(rs)
|
||||
|
||||
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
|
||||
|
|
|
|||
|
|
@ -59,21 +59,21 @@ subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
|
|||
|
||||
ra = max(0d0,rho(iG,1))
|
||||
rb = max(0d0,rho(iG,2))
|
||||
r = ra + rb
|
||||
z = (ra - rb)/r
|
||||
|
||||
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
|
||||
x = sqrt(rs)
|
||||
|
||||
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
|
||||
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
|
||||
|
||||
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
|
||||
|
||||
! spin-up contribution
|
||||
|
||||
if(ra > threshold) then
|
||||
|
||||
r = ra + rb
|
||||
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
|
||||
z = (ra - rb)/r
|
||||
x = sqrt(rs)
|
||||
|
||||
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
|
||||
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
|
||||
|
||||
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
|
||||
|
||||
x_p = x*x + b_p*x + c_p
|
||||
x_f = x*x + b_f*x + c_f
|
||||
x_a = x*x + b_a*x + c_a
|
||||
|
|
@ -132,16 +132,6 @@ subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
|
|||
|
||||
if(rb > threshold) then
|
||||
|
||||
r = ra + rb
|
||||
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
|
||||
z = (ra - rb)/r
|
||||
x = sqrt(rs)
|
||||
|
||||
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
|
||||
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
|
||||
|
||||
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
|
||||
|
||||
x_p = x*x + b_p*x + c_p
|
||||
x_f = x*x + b_f*x + c_f
|
||||
x_a = x*x + b_a*x + c_a
|
||||
|
|
|
|||
|
|
@ -30,6 +30,10 @@ subroutine unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,dr
|
|||
|
||||
call ULYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec)
|
||||
|
||||
case ('PBE')
|
||||
|
||||
call UPBE_gga_correlation_energy(nGrid,weight,rho,drho,Ec)
|
||||
|
||||
case default
|
||||
|
||||
call print_warning('!!! GGA correlation energy not available !!!')
|
||||
|
|
|
|||
|
|
@ -34,7 +34,7 @@ subroutine unrestricted_gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBa
|
|||
|
||||
case ('PBE')
|
||||
|
||||
! call UPBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
|
||||
call UPBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
|
||||
|
||||
case default
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user