LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

Vxc in GW

This commit is contained in:
Pierre-Francois Loos 2021-02-15 17:27:06 +01:00
parent de4927aad4
commit 0f28af6d50
48 changed files with 487 additions and 311 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
include/parameters.h
View File

@ -1,7 +1,7 @@
integer,parameter :: ncart = 3
integer,parameter :: nspin = 2
integer,parameter :: nsp = 3
integer,parameter :: maxEns = 10
integer,parameter :: maxEns = 3
integer,parameter :: maxShell = 512
integer,parameter :: maxL = 7
integer,parameter :: n1eInt = 3
@ -17,7 +17,6 @@
double precision,parameter :: BoToAn = 0.529177249d0
double precision,parameter :: auToD = 2.5415802529d0
double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0)
double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0)
double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0)
double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0)
double precision,parameter :: CxLSDA = - (3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)

6
input/dft
View File

@ -6,18 +6,18 @@
# GGA = 2: B88,G96,PBE
# MGGA = 3:
# Hybrid = 4: HF,B3,PBE
4 HF
4 BHHLYP
# correlation rung:
# Hartree = 0: H
# LDA = 1: VWN5,eVWN5
# GGA = 2: LYP,PBE
# MGGA = 3:
# Hybrid = 4: HF,LYP,PBE
4 HF
4 BHHLYP
# quadrature grid SG-n
1
# Number of states in ensemble (nEns)
3
1
# occupation numbers
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

6
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF KS MOM
F F T F
T F F F
# MP2* MP3 MP2-F12
F F F
# CCD DCD CCSD CCSD(T)
@ -9,11 +9,11 @@
# CIS* CIS(D) CID CISD
F F F F
# RPA* RPAx* ppRPA
F T F
F F F
# G0F2 evGF2 G0F3 evGF3
F F F F
# G0W0* evGW* qsGW*
F F F
T T F
# G0T0 evGT qsGT
F F F
# MCMP2

2
input/options
View File

@ -9,7 +9,7 @@
# GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.0 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
256 0.0000001 T 5 T 0.00367493 F F F F F
256 0.0000001 T 5 T 0.000 F F F F F
# ACFDT: AC Kx XBS
F F T
# BSE: BSE dBSE dTDA evDyn

6
src/AOtoMO/exchange_matrix_AO_basis.f90
View File

@ -1,4 +1,4 @@
subroutine exchange_matrix_AO_basis(nBas,P,G,K)
subroutine exchange_matrix_AO_basis(nBas,P,ERI,K)
! Compute exchange matrix in the AO basis
@ -9,7 +9,7 @@ subroutine exchange_matrix_AO_basis(nBas,P,G,K)
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: G(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
@ -24,7 +24,7 @@ subroutine exchange_matrix_AO_basis(nBas,P,G,K)
do si=1,nBas
do la=1,nBas
do mu=1,nBas
K(mu,nu) = K(mu,nu) - P(la,si)*G(mu,la,si,nu)
K(mu,nu) = K(mu,nu) - P(la,si)*ERI(mu,la,si,nu)
enddo
enddo
enddo

9
src/HF/RHF.f90
View File

@ -1,4 +1,4 @@
subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P)
subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P,Vx)
! Perform restricted Hartree-Fock calculation
@ -55,6 +55,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
double precision,intent(out) :: e(nBas)
double precision,intent(out) :: c(nBas,nBas)
double precision,intent(out) :: P(nBas,nBas)
double precision,intent(out) :: Vx(nBas)
! Hello world
@ -196,7 +197,13 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
ERHF = ET + EV + EJ + EK
! Compute dipole moments
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_RHF(nBas,nO,e,C,ENuc,ET,EV,EJ,EK,ERHF,dipole)
! Compute Vx for post-HF calculations
call exchange_potential(nBas,c,K,Vx)
end subroutine RHF

9
src/HF/UHF.f90
View File

@ -1,4 +1,4 @@
subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P)
subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P,Vx)
! Perform unrestricted Hartree-Fock calculation
@ -59,6 +59,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO
double precision,intent(out) :: e(nBas,nspin)
double precision,intent(out) :: c(nBas,nBas,nspin)
double precision,intent(out) :: P(nBas,nBas,nspin)
double precision,intent(out) :: Vx(nBas,nspin)
! Hello world
@ -249,4 +250,10 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO
call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
! Compute Vx for post-HF calculations
do ispin=1,nspin
call exchange_potential(nBas,c(:,:,ispin),K(:,:,ispin),Vx(:,ispin))
end do
end subroutine UHF

34
src/HF/exchange_potential.f90 Normal file
View File

@ -0,0 +1,34 @@
subroutine exchange_potential(nBas,c,Fx,Vx)
! Compute the exchange potential in the MO basis
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: Fx(nBas,nBas)
! Local variables
integer :: mu,nu
integer :: p
! Output variables
double precision,intent(out) :: Vx(nBas)
! Compute Vx
Vx(:) = 0d0
do p=1,nBas
do mu=1,nBas
do nu=1,nBas
Vx(p) = Vx(p) + c(mu,p)*Fx(mu,nu)*c(nu,p)
end do
end do
end do
end subroutine exchange_potential

24
src/MBPT/G0W0.f90
View File

@ -1,6 +1,6 @@
subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet,triplet,linearize,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,eGW)
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGW)
! Perform G0W0 calculation
@ -29,9 +29,13 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: Vxc(nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
! Local variables
@ -41,6 +45,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcGM
double precision,allocatable :: SigX(:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: OmRPA(:)
@ -100,14 +105,14 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
! Memory allocation
allocate(SigC(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGWlin(nBas))
allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGWlin(nBas))
!-------------------!
! Compute screening !
!-------------------!
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0, &
eHF,ERI,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
eHF,ERI_MO,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA)
@ -115,12 +120,13 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
! Compute spectral weights !
!--------------------------!
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
!------------------------!
! Compute GW self-energy !
!------------------------!
call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX)
call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
!--------------------------------!
@ -133,7 +139,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
! Solve the quasi-particle equation !
!-----------------------------------!
eGWlin(:) = eHF(:) + Z(:)*SigC(:)
eGWlin(:) = eHF(:) + Z(:)*(SigX(:) + SigC(:) - Vxc(:))
! Linearized or graphical solution?
@ -159,7 +165,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
! Compute the RPA correlation energy
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, &
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
!--------------!
@ -180,7 +186,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
if(BSE) then
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE)
if(exchange_kernel) then
@ -214,7 +220,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcAC)
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

25
src/MBPT/UG0W0.f90
View File

@ -1,6 +1,6 @@
subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eGW)
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eGW)
! Perform unrestricted G0W0 calculation
@ -34,6 +34,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
@ -42,6 +43,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: PHF(nBas,nBas,nspin)
double precision,intent(in) :: Vxc(nBas,nspin)
! Local variables
@ -52,6 +54,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
double precision :: EcGM(nspin)
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision,allocatable :: SigX(:,:)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: Z(:,:)
integer :: nS_aa,nS_bb,nS_sc
@ -106,8 +109,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(SigC(nBas,nspin),Z(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc), &
rho_RPA(nBas,nBas,nS_sc,nspin),eGWlin(nBas,nspin))
allocate(SigX(nBas,nspin),SigC(nBas,nspin),Z(nBas,nspin),eGWlin(nBas,nspin), &
OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin))
!-------------------!
! Compute screening !
@ -132,6 +135,10 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! Compute self-energy !
!---------------------!
do is=1,nspin
call self_energy_exchange_diag(nBas,cHF(:,:,is),PHF(:,:,is),ERI,SigX(:,is))
end do
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM)
!--------------------------------!
@ -144,7 +151,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! Solve the quasi-particle equation !
!-----------------------------------!
eGWlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:)
eGWlin(:,:) = eHF(:,:) + Z(:,:)*(SigX(:,:) + SigC(:,:) - Vxc(:,:))
if(linearize) then
@ -157,10 +164,10 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! Find graphical solution of the QP equation
do is=1,nspin
call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),OmRPA, &
rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is))
end do
do is=1,nspin
call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),OmRPA, &
rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is))
end do
end if

25
src/MBPT/evGW.f90
View File

@ -1,6 +1,6 @@
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI, &
dipole_int,eHF,eG0W0)
G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0)
! Perform self-consistent eigenvalue-only GW calculation
@ -31,9 +31,13 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: Vxc(nBas,nspin)
double precision,intent(in) :: eG0W0(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
@ -54,6 +58,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
double precision,allocatable :: eGW(:)
double precision,allocatable :: eOld(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: SigX(:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
@ -117,9 +122,13 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
! Memory allocation
allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), &
allocate(eGW(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), &
rho_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Compute the exchange part of the self-energy
call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX)
! Initialization
nSCF = 0
@ -142,14 +151,14 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
if(.not. GW0 .or. nSCF == 0) then
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, &
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
endif
! Compute spectral weights
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
! Compute correlation part of the self-energy
@ -233,7 +242,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
if(BSE) then
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE)
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eGW,eGW,EcBSE)
if(exchange_kernel) then
@ -267,7 +276,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcAC)
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

19
src/MBPT/evUGW.f90
View File

@ -1,6 +1,6 @@
subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, &
EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eG0W0)
EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0)
! Perform self-consistent eigenvalue-only GW calculation
@ -37,10 +37,13 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: PHF(nBas,nBas,nspin)
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: Vxc(nBas,nspin)
double precision,intent(in) :: eG0W0(nBas,nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
@ -67,6 +70,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
double precision,allocatable :: eOld(:,:)
double precision,allocatable :: Z(:,:)
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: SigX(:,:)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
@ -127,8 +131,15 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigC(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc), &
XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin),error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin))
allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigX(nBas,nspin),SigC(nBas,nspin), &
OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin), &
error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin))
! Compute the exchange part of the self-energy
do is=1,nspin
call self_energy_exchange_diag(nBas,cHF(:,:,is),PHF(:,:,is),ERI_AO,SigX(:,is))
end do
! Initialization
@ -182,7 +193,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
! Solve the quasi-particle equation !
!-----------------------------------!
eGW(:,:) = eHF(:,:) + SigC(:,:)
eGW(:,:) = eHF(:,:) + SigX(:,:) + SigC(:,:) - Vxc(:,:)
! Convergence criteria

10
src/MBPT/self_energy_exchange.f90
View File

@ -1,4 +1,4 @@
subroutine self_energy_exchange(nBas,c,P,G,SigmaX)
subroutine self_energy_exchange(nBas,c,P,ERI,SigX)
! Compute exchange part of the self-energy
@ -8,18 +8,18 @@ subroutine self_energy_exchange(nBas,c,P,G,SigmaX)
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas),G(nBas,nBas,nBas,nBas)
double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas),ERI(nBas,nBas,nBas,nBas)
! Output variables
double precision,intent(out) :: SigmaX(nBas,nBas)
double precision,intent(out) :: SigX(nBas,nBas)
! Compute exchange part of the self-energy in the AO basis
call exchange_matrix_AO_basis(nBas,P,G,SigmaX)
call exchange_matrix_AO_basis(nBas,P,ERI,SigX)
! Compute exchange part of the self-energy in the MO basis
SigmaX = matmul(transpose(c),matmul(SigmaX,c))
SigX = matmul(transpose(c),matmul(SigX,c))
end subroutine self_energy_exchange

42
src/MBPT/self_energy_exchange_diag.f90 Normal file
View File

@ -0,0 +1,42 @@
subroutine self_energy_exchange_diag(nBas,c,P,ERI,SigX)
! Compute the diagonal elements of the exchange part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: q,mu,nu
double precision,allocatable :: Fx(:,:)
! Output variables
double precision,intent(out) :: SigX(nBas)
! Compute exchange part of the self-energy in the AO basis
allocate(Fx(nBas,nBas))
call exchange_matrix_AO_basis(nBas,P,ERI,Fx)
! Compute exchange part of the self-energy in the MO basis
SigX(:) = 0d0
do q=1,nBas
do mu=1,nBas
do nu=1,nBas
SigX(q) = SigX(q) + c(mu,q)*Fx(mu,nu)*c(nu,q)
end do
end do
end do
deallocate(Fx)
end subroutine self_energy_exchange_diag

104
src/QuAcK/QuAcK.f90
View File

@ -32,6 +32,7 @@ program QuAcK
double precision,allocatable :: ZNuc(:),rNuc(:,:)
double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:)
double precision,allocatable :: Vxc(:,:,:)
double precision,allocatable :: eG0W0(:,:)
double precision,allocatable :: eG0T0(:,:)
@ -236,7 +237,7 @@ program QuAcK
allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), &
S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas), &
dipole_int_AO(nBas,nBas,ncart),dipole_int_MO(nBas,nBas,ncart))
dipole_int_AO(nBas,nBas,ncart),dipole_int_MO(nBas,nBas,ncart),Vxc(nBas,nBas,nspin))
! Read integrals
@ -280,7 +281,7 @@ program QuAcK
call cpu_time(start_HF)
call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF)
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc)
call cpu_time(end_HF)
t_HF = end_HF - start_HF
@ -300,7 +301,7 @@ program QuAcK
call cpu_time(start_HF)
call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF)
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF,Vxc)
call cpu_time(end_HF)
t_HF = end_HF - start_HF
@ -315,10 +316,13 @@ program QuAcK
if(doKS) then
! Switch on the unrestricted flag
unrestricted = .true.
call cpu_time(start_KS)
call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, &
nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EUHF,eHF,cHF,PHF)
max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EUHF,eHF,cHF,PHF,Vxc)
call cpu_time(end_KS)
@ -851,12 +855,12 @@ program QuAcK
if(unrestricted) then
call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0)
linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0)
else
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,eG0W0)
linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
end if
@ -879,14 +883,14 @@ program QuAcK
call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc, &
EUHF,S,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, &
cHF,eHF,eG0W0)
EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, &
PHF,cHF,eHF,Vxc,eG0W0)
else
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, &
BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,eG0W0)
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
end if
call cpu_time(end_evGW)
@ -1005,10 +1009,10 @@ program QuAcK
if(doMCMP2) then
call cpu_time(start_MCMP2)
call MCMP2(doDrift,nBas,nC,nO,nV,cHF,eHF,EcMP2, &
nMC,nEq,nWalk,dt,nPrint, &
nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
Norm,EcMCMP2,Err_EcMCMP2,Var_EcMCMP2)
! call MCMP2(doDrift,nBas,nC,nO,nV,cHF,eHF,EcMP2, &
! nMC,nEq,nWalk,dt,nPrint, &
! nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
! Norm,EcMCMP2,Err_EcMCMP2,Var_EcMCMP2)
call cpu_time(end_MCMP2)
t_MCMP2 = end_MCMP2 - start_MCMP2
@ -1040,59 +1044,59 @@ program QuAcK
! Range-separeted GT/GW
!------------------------------------------------------------------------
if(doGTGW) then
! if(doGTGW) then
! Read and transform long-range two-electron integrals
! ! Read and transform long-range two-electron integrals
allocate(ERI_ERF_AO(nBas,nBas,nBas,nBas),ERI_ERF_MO(nBas,nBas,nBas,nBas))
call read_LR(nBas,ERI_ERF_AO)
! allocate(ERI_ERF_AO(nBas,nBas,nBas,nBas),ERI_ERF_MO(nBas,nBas,nBas,nBas))
! call read_LR(nBas,ERI_ERF_AO)
call cpu_time(start_AOtoMO)
! call cpu_time(start_AOtoMO)
write(*,*)
write(*,*) 'AO to MO transformation for long-range ERIs... Please be patient'
write(*,*)
! write(*,*)
! write(*,*) 'AO to MO transformation for long-range ERIs... Please be patient'
! write(*,*)
call AOtoMO_integral_transform(nBas,cHF,ERI_ERF_AO,ERI_ERF_MO)
! call AOtoMO_integral_transform(nBas,cHF,ERI_ERF_AO,ERI_ERF_MO)
call cpu_time(end_AOtoMO)
! call cpu_time(end_AOtoMO)
deallocate(ERI_ERF_AO)
! deallocate(ERI_ERF_AO)
t_AOtoMO = end_AOtoMO - start_AOtoMO
! t_AOtoMO = end_AOtoMO - start_AOtoMO
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds'
write(*,*)
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds'
! write(*,*)
! Long-range G0W0 calculation
! ! Long-range G0W0 calculation
call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0W0)
call cpu_time(end_G0W0)
t_G0W0 = end_G0W0 - start_G0W0
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_G0W0,' seconds'
write(*,*)
! call cpu_time(start_G0W0)
! call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
! dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, &
! nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0W0)
! call cpu_time(end_G0W0)
!
! t_G0W0 = end_G0W0 - start_G0W0
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_G0W0,' seconds'
! write(*,*)
! Short-range G0T0 calculation
! ! Short-range G0T0 calculation
ERI_ERF_MO(:,:,:,:) = ERI_MO(:,:,:,:) - ERI_ERF_MO(:,:,:,:)
! ERI_ERF_MO(:,:,:,:) = ERI_MO(:,:,:,:) - ERI_ERF_MO(:,:,:,:)
call cpu_time(start_G0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, &
singlet,triplet,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0T0)
call cpu_time(end_G0T0)
t_G0T0 = end_G0T0 - start_G0T0
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds'
write(*,*)
! call cpu_time(start_G0T0)
! call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, &
! singlet,triplet,linGW,eta_GW, &
! nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0T0)
! call cpu_time(end_G0T0)
!
! t_G0T0 = end_G0T0 - start_G0T0
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds'
! write(*,*)
! call matout(nBas,1,(eG0W0+eG0T0-eHF(:,1))*HaToeV)
end if
! end if
!------------------------------------------------------------------------
! Basis set correction

4
src/eDFT/RB88_gga_exchange_energy.f90
View File

@ -16,7 +16,6 @@ subroutine RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Local variables
integer :: iG
double precision :: alpha
double precision :: beta
double precision :: r,g,x
@ -26,7 +25,6 @@ subroutine RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Coefficients for B88 GGA exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
beta = 0.0042d0
! Compute GGA exchange energy
@ -41,7 +39,7 @@ subroutine RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
g = 0.25d0*(drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2)
x = sqrt(g)/r**(4d0/3d0)
Ex = Ex + weight(iG)*alpha*r**(4d0/3d0) &
Ex = Ex + weight(iG)*CxLDA*r**(4d0/3d0) &
- weight(iG)*beta*x**2*r**(4d0/3d0)/(1d0 + 6d0*beta*x*asinh(x))
end if

6
src/eDFT/RB88_gga_exchange_individual_energy.f90
View File

@ -17,7 +17,6 @@ subroutine RB88_gga_exchange_individual_energy(nGrid,weight,rhow,drhow,rho,drho,
! Local variables
integer :: iG
double precision :: alpha
double precision :: beta
double precision :: r,rI,g,x
double precision :: ex_p,dexdr_p
@ -28,7 +27,6 @@ subroutine RB88_gga_exchange_individual_energy(nGrid,weight,rhow,drhow,rho,drho,
! Coefficients for B88 GGA exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
beta = 0.0042d0
! Compute GGA exchange matrix in the AO basis
@ -45,11 +43,11 @@ subroutine RB88_gga_exchange_individual_energy(nGrid,weight,rhow,drhow,rho,drho,
g = 0.25d0*(drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2)
x = sqrt(g)/r**(4d0/3d0)
dexdr_p = 4d0/3d0*r**(1d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2) &
dexdr_p = 4d0/3d0*r**(1d0/3d0)*(CxLDA - beta*g**(3d0/4d0)/r**2) &
+ 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0) &
- 2d0*3d0/4d0*beta*g**(-1d0/4d0)/r**(2d0/3d0)
ex_p = alpha*r**(4d0/3d0) &
ex_p = CxLDA*r**(4d0/3d0) &
- weight(iG)*beta*x**2*r**(4d0/3d0)/(1d0 + 6d0*beta*x*asinh(x))
Ex = Ex + weight(iG)*(ex_p*rI + dexdr_p*r*rI - dexdr_p*r*r)

4
src/eDFT/RB88_gga_exchange_potential.f90
View File

@ -18,7 +18,6 @@ subroutine RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Local variables
integer :: mu,nu,iG
double precision :: alpha
double precision :: beta
double precision :: r,g,vAO,gAO
@ -28,7 +27,6 @@ subroutine RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Coefficients for B88 GGA exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
beta = 0.0042d0
! Compute GGA exchange matrix in the AO basis
@ -46,7 +44,7 @@ subroutine RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
g = 0.25d0*(drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2)
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) &
+ vAO*(4d0/3d0*r**(1d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2) &
+ vAO*(4d0/3d0*r**(1d0/3d0)*(CxLDA - beta*g**(3d0/4d0)/r**2) &
+ 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0))
gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &

6
src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90
View File

@ -84,12 +84,16 @@ subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)) &
* (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
case default
dCxdw1 = 0.d0
dCxdw2 = 0.d0
end select
dCxdw1 = CxLDA*dCxdw1
dCxdw2 = CxLDA*dCxdw2
dExdw(:) = 0d0
do iG=1,nGrid

19
src/eDFT/RCC_lda_exchange_energy.f90
View File

@ -69,12 +69,19 @@ subroutine RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,C
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
select case (Cx_choice)
case(1)
Cx = CxLDA*Fx1
case(2)
Cx = CxLDA*Fx2
case(3)
Cx = CxLDA*Fx2*Fx1
case(1)
Cx = CxLDA*Fx1
case(2)
Cx = CxLDA*Fx2
case(3)
Cx = CxLDA*Fx2*Fx1
case default
Cx = CxLDA
end select
! Compute GIC-LDA exchange energy

19
src/eDFT/RCC_lda_exchange_individual_energy.f90
View File

@ -75,12 +75,19 @@ subroutine RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
select case (Cx_choice)
case(1)
Cx = CxLDA*Fx1
case(2)
Cx = CxLDA*Fx2
case(3)
Cx = CxLDA*Fx2*Fx1
case(1)
Cx = CxLDA*Fx1
case(2)
Cx = CxLDA*Fx2
case(3)
Cx = CxLDA*Fx2*Fx1
case default
Cx = CxLDA
end select
! Compute LDA exchange matrix in the AO basis

19
src/eDFT/RCC_lda_exchange_potential.f90
View File

@ -74,12 +74,19 @@ subroutine RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
select case (Cx_choice)
case(1)
Cx = CxLDA*Fx1
case(2)
Cx = CxLDA*Fx2
case(3)
Cx = CxLDA*Fx2*Fx1
case(1)
Cx = CxLDA*Fx1
case(2)
Cx = CxLDA*Fx2
case(3)
Cx = CxLDA*Fx2*Fx1
case default
Cx = CxLDA
end select

3
src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90
View File

@ -21,6 +21,9 @@ subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,r
double precision :: dExdw(nEns)
double precision,external :: Kronecker_delta
double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0)
double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0)
! Output variables
double precision,intent(out) :: ExDD(nEns)

3
src/eDFT/RMFL20_lda_exchange_energy.f90
View File

@ -15,6 +15,9 @@ subroutine RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0)
double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0)
! Local variables
integer :: iG

3
src/eDFT/RMFL20_lda_exchange_individual_energy.f90
View File

@ -22,6 +22,9 @@ subroutine RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,we
double precision :: r,rI
double precision :: e_p,dedr
double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0)
double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0)
! Output variables
double precision,intent(out) :: Ex

3
src/eDFT/RMFL20_lda_exchange_potential.f90
View File

@ -22,6 +22,9 @@ subroutine RMFL20_lda_exchange_potential(LDA_centered,nEns,wEns,nGrid,weight,nBa
double precision :: Cxw
double precision :: r,vAO
double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0)
double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0)
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)

5
src/eDFT/UB88_gga_exchange_energy.f90
View File

@ -16,7 +16,7 @@ subroutine UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Local variables
integer :: iG
double precision :: alpha,b
double precision :: b
double precision :: r,g,x
! Output variables
@ -25,7 +25,6 @@ subroutine UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Coefficients for B88 GGA exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
b = 0.0042d0
! Compute GGA exchange energy
@ -40,7 +39,7 @@ subroutine UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
x = sqrt(g)/r**(4d0/3d0)
Ex = Ex + weight(iG)*r**(4d0/3d0)*(alpha - b*x**2/(1d0 + 6d0*b*x*asinh(x)))
Ex = Ex + weight(iG)*r**(4d0/3d0)*(CxLSDA - b*x**2/(1d0 + 6d0*b*x*asinh(x)))
end if

9
src/eDFT/UB88_gga_exchange_potential.f90
View File

@ -18,7 +18,7 @@ subroutine UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Local variables
integer :: mu,nu,iG
double precision :: alpha,b
double precision :: b
double precision :: vAO,gAO
double precision :: r,g,x,dxdr,dxdg,f
@ -28,7 +28,6 @@ subroutine UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Coefficients for B88 GGA exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
b = 0.0042d0
! Compute GGA exchange matrix in the AO basis
@ -52,9 +51,9 @@ subroutine UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
f = b*x**2/(1d0 + 6d0*b*x*asinh(x))
Fx(mu,nu) = Fx(mu,nu) + vAO*( &
4d0/3d0*r**(1d0/3d0)*(alpha - f) &
- 2d0*r**(4d0/3d0)*dxdr*f &
Fx(mu,nu) = Fx(mu,nu) + vAO*( &
4d0/3d0*r**(1d0/3d0)*(CxLSDA - f) &
- 2d0*r**(4d0/3d0)*dxdr*f &
+ r**(4d0/3d0)*dxdr*(6d0*b*x*asinh(x) + 6d0*b*x**2/sqrt(1d0+x**2))*f/(1d0 + 6d0*b*x*asinh(x)) )
gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &

22
src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90
View File

@ -22,7 +22,7 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
integer :: iEns,jEns
integer :: iG
double precision :: r,alpha
double precision :: r
double precision,allocatable :: dExdw(:)
double precision,external :: Kronecker_delta
@ -30,9 +30,6 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
double precision :: a2,b2,c2,w2
double precision :: dCxdw1,dCxdw2
double precision :: nEli,nElw
! Output variables
double precision,intent(out) :: ExDD(nEns)
@ -60,13 +57,6 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
b2 = aCC_w2(2)
c2 = aCC_w2(3)
nElw = electron_number(nGrid,weight,rhow)
! Cx coefficient for unrestricted Slater LDA exchange
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
w1 = wEns(2)
w2 = wEns(3)
@ -86,11 +76,15 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)) &
* (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
case default
dCxdw1 = 0d0
dCxdw2 = 0d0
end select
dCxdw1 = alpha*dCxdw1
dCxdw2 = alpha*dCxdw2
dCxdw1 = CxLSDA*dCxdw1
dCxdw2 = CxLSDA*dCxdw2
dExdw(:) = 0d0

33
src/eDFT/UCC_lda_exchange_energy.f90
View File

@ -19,7 +19,7 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,C
! Local variables
integer :: iG
double precision :: r,alpha
double precision :: r
double precision :: a1,b1,c1,w1
double precision :: a2,b2,c2,w2
@ -65,10 +65,6 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,C
b2 = aCC_w2(2)
c2 = aCC_w2(3)
! Cx coefficient for unrestricted Slater LDA exchange
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
! Fx1 for states N and N-1
! Fx2 for states N and N+1
@ -79,21 +75,20 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,C
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
select case (Cx_choice)
case(1)
Cx = alpha*Fx1
case(2)
Cx = alpha*Fx2
case(3)
Cx = alpha*Fx2*Fx1
case(1)
Cx = CxLSDA*Fx1
case(2)
Cx = CxLSDA*Fx2
case(3)
Cx = CxLSDA*Fx2*Fx1
case default
Cx = CxLSDA
end select
! for two-weights ensemble
! Cx = alpha*Fx2*Fx1
! for left ensemble
! Cx = alpha*Fx1
! for right ensemble
! Cx = alpha*Fx2
! Compute GIC-LDA exchange energy

34
src/eDFT/UCC_lda_exchange_individual_energy.f90
View File

@ -22,9 +22,8 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
! Local variables
integer :: iG
double precision :: r,rI,alpha
double precision :: r,rI
double precision :: e_p,dedr
double precision :: nEli,nElw
double precision :: a1,b1,c1,w1
double precision :: a2,b2,c2,w2
@ -51,10 +50,6 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
b2 = aCC_w2(2)
c2 = aCC_w2(3)
! Cx coefficient for unrestricted Slater LDA exchange
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
w1 = wEns(2)
Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)
@ -64,21 +59,19 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
select case (Cx_choice)
case(1)
Cx = alpha*Fx1
Cx = CxLSDA*Fx1
case(2)
Cx = alpha*Fx2
Cx = CxLSDA*Fx2
case(3)
Cx = alpha*Fx2*Fx1
Cx = CxLSDA*Fx2*Fx1
case default
Cx = CxLSDA
end select
nEli = electron_number(nGrid,weight,rho)
nElw = electron_number(nGrid,weight,rhow)
! Compute LDA exchange matrix in the AO basis
Ex = 0d0
@ -92,19 +85,12 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
e_p = Cx*r**(1d0/3d0)
dedr = 1d0/3d0*Cx*r**(-2d0/3d0)
if (doNcentered) then
Ex = Ex - weight(iG)*dedr*r*r*(nEli/nElw)
else
Ex = Ex - weight(iG)*dedr*r*r
end if
Ex = Ex - weight(iG)*dedr*r*r
if(rI > threshold) then
if (doNcentered) then
Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI)
else
Ex = Ex + weight(iG)*((nEli/nElw)*e_p*rI + dedr*r*rI)
end if
Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI)
endif
endif

34
src/eDFT/UCC_lda_exchange_potential.f90
View File

@ -21,7 +21,7 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
! Local variables
integer :: mu,nu,iG
double precision :: r,vAO,alpha
double precision :: r,vAO
double precision :: a1,b1,c1,w1
double precision :: a2,b2,c2,w2
@ -67,10 +67,6 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
b2 = aCC_w2(2)
c2 = aCC_w2(3)
! Cx coefficient for unrestricted Slater LDA exchange
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
! Fx1 for states N and N-1
! Fx2 for states N and N+1
@ -81,23 +77,21 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
select case (Cx_choice)
case(1)
Cx = alpha*Fx1
case(2)
Cx = alpha*Fx2
case(3)
Cx = alpha*Fx2*Fx1
case(1)
Cx = CxLSDA*Fx1
case(2)
Cx = CxLSDA*Fx2
case(3)
Cx = CxLSDA*Fx2*Fx1
case default
Cx = CxLSDA
end select
! for two-weight ensembles
! Cx = alpha*Fx2*Fx1
! for left ensemble
! Cx = alpha*Fx1
! for right ensemble
! Cx = alpha*Fx2
! Compute LDA exchange matrix in the AO basis
Fx(:,:) = 0d0

5
src/eDFT/UG96_gga_exchange_energy.f90
View File

@ -16,7 +16,7 @@ subroutine UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Local variables
integer :: iG
double precision :: alpha,beta
double precision :: beta
double precision :: r,g
! Output variables
@ -25,7 +25,6 @@ subroutine UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Coefficients for G96 GGA exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
beta = 1d0/137d0
! Compute GGA exchange energy
@ -40,7 +39,7 @@ subroutine UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
Ex = Ex + weight(iG)*r**(4d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2)
Ex = Ex + weight(iG)*r**(4d0/3d0)*(CxLSDA - beta*g**(3d0/4d0)/r**2)
end if

5
src/eDFT/UG96_gga_exchange_potential.f90
View File

@ -18,7 +18,7 @@ subroutine UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Local variables
integer :: mu,nu,iG
double precision :: alpha,beta
double precision :: beta
double precision :: r,g,vAO,gAO
! Output variables
@ -27,7 +27,6 @@ subroutine UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Coefficients for G96 GGA exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
beta = +1d0/137d0
! Compute GGA exchange matrix in the AO basis
@ -45,7 +44,7 @@ subroutine UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) &
+ vAO*(4d0/3d0*r**(1d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2) &
+ vAO*(4d0/3d0*r**(1d0/3d0)*(CxLSDA - beta*g**(3d0/4d0)/r**2) &
+ 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0))
gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &

5
src/eDFT/UPBE_gga_exchange_energy.f90
View File

@ -16,7 +16,7 @@ subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Local variables
integer :: iG
double precision :: alpha,mupbe,kappa
double precision :: mupbe,kappa
double precision :: r,g,s2
! Output variables
@ -25,7 +25,6 @@ subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Coefficients for PBE exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
mupbe = ((1d0/2d0**(1d0/3d0))/(2d0*(3d0*pi**2)**(1d0/3d0)))**2*0.21951d0
kappa = 0.804d0
@ -41,7 +40,7 @@ subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
s2 = g/r**(8d0/3d0)
Ex = Ex + weight(iG)*alpha*r**(4d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa))
Ex = Ex + weight(iG)*CxLSDA*r**(4d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa))
end if

9
src/eDFT/UPBE_gga_exchange_potential.f90
View File

@ -18,7 +18,7 @@ subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Local variables
integer :: mu,nu,iG
double precision :: alpha,mupbe,kappa
double precision :: mupbe,kappa
double precision :: r,g,s2,vAO,gAO
! Output variables
@ -27,7 +27,6 @@ subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Coefficients for PBE exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
mupbe = ((1d0/2d0**(1d0/3d0))/(2d0*(3d0*pi**2)**(1d0/3d0)))**2*0.21951d0
kappa = 0.804d0
@ -49,15 +48,15 @@ subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) &
+ vAO*4d0/3d0*alpha*r**(1d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa)) &
- vAO*8d0/3d0*alpha*r**(1d0/3d0)*mupbe*s2/(1d0 + mupbe*s2/kappa)**2
+ vAO*4d0/3d0*CxLSDA*r**(1d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa)) &
- vAO*8d0/3d0*CxLSDA*r**(1d0/3d0)*mupbe*s2/(1d0 + mupbe*s2/kappa)**2
gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &
+ drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) &
+ drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG))
gAO = weight(iG)*gAO
Fx(mu,nu) = Fx(mu,nu) + 2d0*gAO*alpha*r**(-4d0/3d0)*mupbe/(1d0 + mupbe*s2/kappa)**2
Fx(mu,nu) = Fx(mu,nu) + 2d0*gAO*CxLSDA*r**(-4d0/3d0)*mupbe/(1d0 + mupbe*s2/kappa)**2
end if

8
src/eDFT/US51_lda_exchange_energy.f90
View File

@ -16,16 +16,12 @@ subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex)
! Local variables
integer :: iG
double precision :: alpha,r,alphaw,a2,b2,c2,a1,b1,c1
double precision :: r
! Output variables
double precision :: Ex
! Cx coefficient for Slater LDA exchange
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
! Compute LDA exchange energy
Ex = 0d0
@ -35,7 +31,7 @@ subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex)
if(r > threshold) then
Ex = Ex + weight(iG)*alpha*r**(4d0/3d0)
Ex = Ex + weight(iG)*CxLSDA*r**(4d0/3d0)
endif

8
src/eDFT/US51_lda_exchange_individual_energy.f90
View File

@ -17,7 +17,7 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered
! Local variables
integer :: iG
double precision :: r,rI,alpha
double precision :: r,rI
double precision :: e,dedr
double precision :: Exrr,ExrI,ExrrI
@ -27,8 +27,6 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered
! Compute LDA exchange matrix in the AO basis
alpha = - (3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
Exrr = 0d0
ExrI = 0d0
ExrrI = 0d0
@ -40,8 +38,8 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered
if(r > threshold) then
e = alpha*r**(1d0/3d0)
dedr = 1d0/3d0*alpha*r**(-2d0/3d0)
e = CxLSDA*r**(1d0/3d0)
dedr = 1d0/3d0*CxLSDA*r**(-2d0/3d0)
Exrr = Exrr - weight(iG)*dedr*r*r

7
src/eDFT/US51_lda_exchange_potential.f90
View File

@ -16,17 +16,12 @@ subroutine US51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
! Local variables
integer :: mu,nu,iG
double precision :: alpha
double precision :: r,vAO
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Cx coefficient for Slater LDA exchange
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
! Compute LDA exchange matrix in the AO basis
Fx(:,:) = 0d0
@ -39,7 +34,7 @@ subroutine US51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
if(r > threshold) then
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*alpha*r**(1d0/3d0)
Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*CxLSDA*r**(1d0/3d0)
endif

4
src/eDFT/UVWN5_lda_correlation_individual_energy.f90
View File

@ -88,7 +88,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
decdra_f = drsdra*dxdrs*decdx_f
Ec(1) = Ec(1) + weight(iG)*(ec_z + decdra_f*r)*rI
Ec(1) = Ec(1) + weight(iG)*(ec_f + decdra_f*r)*rI
end if
@ -191,7 +191,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
decdra_f = drsdra*dxdrs*decdx_f
Ec(3) = Ec(3) + weight(iG)*(ec_z + decdra_f*r)*rI
Ec(3) = Ec(3) + weight(iG)*(ec_f + decdra_f*r)*rI
end if

89
src/eDFT/eDFT.f90
View File

@ -1,6 +1,6 @@
subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, &
nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI,dipole_int,Ew,eKS,cKS,PKS)
max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI,dipole_int,Ew,eKS,cKS,PKS,Vxc)
! exchange-correlation density-functional theory calculations
@ -50,8 +50,6 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
! Local variables
double precision,allocatable :: c(:,:)
character(len=8) :: method
integer :: x_rung,c_rung
character(len=12) :: x_DFA ,c_DFA
@ -88,6 +86,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
double precision,intent(out) :: eKS(nBas,nspin)
double precision,intent(out) :: cKS(nBas,nBas,nspin)
double precision,intent(out) :: PKS(nBas,nBas,nspin)
double precision,intent(out) :: Vxc(nBas,nspin)
! Hello World
@ -111,7 +110,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
! Allocate ensemble weights and MO coefficients
allocate(c(nBas,nspin),wEns(maxEns),occnum(nBas,nspin,maxEns))
allocate(wEns(maxEns),occnum(nBas,nspin,maxEns))
call read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, &
doNcentered,occnum,Cx_choice)
@ -140,73 +139,73 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
! Compute GOK-RKS energy
!------------------------------------------------------------------------
if(method == 'GOK-RKS') then
! if(method == 'GOK-RKS') then
call cpu_time(start_KS)
call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, &
maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO(1),nV(1), &
S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum,Cx_choice)
call cpu_time(end_KS)
! call cpu_time(start_KS)
! call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, &
! maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO(1),nV(1), &
! S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum,Cx_choice)
! call cpu_time(end_KS)
t_KS = end_KS - start_KS
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GOK-RKS = ',t_KS,' seconds'
write(*,*)
! t_KS = end_KS - start_KS
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GOK-RKS = ',t_KS,' seconds'
! write(*,*)
end if
!end if
!------------------------------------------------------------------------
! Compute LIM excitation energies
!------------------------------------------------------------------------
if(method == 'LIM-RKS') then
! if(method == 'LIM-RKS') then
call cpu_time(start_KS)
call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), &
aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), &
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered)
call cpu_time(end_KS)
! call cpu_time(start_KS)
! call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), &
! aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), &
! S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered)
! call cpu_time(end_KS)
t_KS = end_KS - start_KS
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for LIM-RKS = ',t_KS,' seconds'
write(*,*)
! t_KS = end_KS - start_KS
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for LIM-RKS = ',t_KS,' seconds'
! write(*,*)
end if
! end if
!------------------------------------------------------------------------
! Compute MOM excitation energies
!------------------------------------------------------------------------
if(method == 'MOM-RKS') then
! if(method == 'MOM-RKS') then
call cpu_time(start_KS)
call MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), &
aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), &
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered)
call cpu_time(end_KS)
! call cpu_time(start_KS)
! call MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), &
! aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), &
! S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered)
! call cpu_time(end_KS)
t_KS = end_KS - start_KS
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM-RKS = ',t_KS,' seconds'
write(*,*)
! t_KS = end_KS - start_KS
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM-RKS = ',t_KS,' seconds'
! write(*,*)
end if
! end if
!------------------------------------------------------------------------
! Compute GOK-UKS energy (BROKEN)
!------------------------------------------------------------------------
if(method == 'GOK-UKS') then
! if(method == 'GOK-UKS') then
call cpu_time(start_KS)
call GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, &
nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum, &
Cx_choice,doNcentered)
call cpu_time(end_KS)
! call cpu_time(start_KS)
! call GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, &
! nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum, &
! Cx_choice,doNcentered)
! call cpu_time(end_KS)
t_KS = end_KS - start_KS
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UKS = ',t_KS,' seconds'
write(*,*)
! t_KS = end_KS - start_KS
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UKS = ',t_KS,' seconds'
! write(*,*)
end if
! end if
!------------------------------------------------------------------------
! Compute N-centered UKS energy
@ -216,7 +215,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
call cpu_time(start_KS)
call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),maxSCF,thresh,max_diis,guess_type,mix, &
nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS)
nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
call cpu_time(end_KS)
t_KS = end_KS - start_KS

16
src/eDFT/eDFT_UKS.f90
View File

@ -1,5 +1,5 @@
subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,occnum,Cx_choice,doNcentered,Ew,eps,c,Pw)
nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,occnum,Cx_choice,doNcentered,Ew,eps,c,Pw,Vxc)
! Perform unrestricted Kohn-Sham calculation for ensembles
@ -80,7 +80,8 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
double precision,intent(out) :: Ew
double precision,intent(out) :: eps(nBas,nspin)
double precision,intent(out) :: Pw(nBas,nBas,nspin)
double precision,intent(out) :: c(nBAs,nBas,nspin)
double precision,intent(out) :: c(nBas,nBas,nspin)
double precision,intent(out) :: Vxc(nBas,nspin)
! Hello world
@ -386,12 +387,17 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
call print_UKS(nBas,nEns,occnum,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew)
! Compute Vxc for post-HF calculations
call xc_potential(nBas,c,Fx,Fc,Vxc)
!------------------------------------------------------------------------
! Compute individual energies from ensemble energy
!------------------------------------------------------------------------
call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
AO,dAO,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum, &
Cx_choice,doNcentered)
if(nEns > 1) &
call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
AO,dAO,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum, &
Cx_choice,doNcentered)
end subroutine eDFT_UKS

6
src/eDFT/read_options_dft.f90
View File

@ -96,7 +96,7 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
occnum(:,:,:) = 0d0
do iEns=1,nEns
do iEns=1,maxEns
read(1,*)
read(1,*) (occnum(iBas,1,iEns),iBas=1,nBas)
read(1,*) (occnum(iBas,2,iEns),iBas=1,nBas)
@ -117,9 +117,9 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
! Read ensemble weights for real physical (fractional number of electrons) ensemble (w1,w2)
allocate(nEl(nEns))
allocate(nEl(maxEns))
nEl(:) = 0d0
do iEns=1,nEns
do iEns=1,maxEns
do iBas=1,nBas
nEl(iEns) = nEl(iEns) + occnum(iBas,1,iEns) + occnum(iBas,2,iEns)
end do

10
src/eDFT/unrestricted_hybrid_correlation_energy.f90
View File

@ -31,7 +31,7 @@ subroutine unrestricted_hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho
Ec(:) = 0d0
case('LYP')
case('B3LYP')
aC = 0.81d0
@ -40,11 +40,13 @@ subroutine unrestricted_hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho
Ec(:) = EcLDA(:) + aC*(EcGGA(:) - EcLDA(:))
case('BHHLYP')
call unrestricted_gga_correlation_energy('LYP ',nEns,wEns,nGrid,weight,rho,drho,Ec)
case('PBE')
call unrestricted_gga_correlation_energy('PBE ',nEns,wEns,nGrid,weight,rho,drho,EcGGA)
Ec(:) = EcGGA(:)
call unrestricted_gga_correlation_energy('PBE ',nEns,wEns,nGrid,weight,rho,drho,Ec)
case default

12
src/eDFT/unrestricted_hybrid_correlation_potential.f90
View File

@ -36,7 +36,7 @@ subroutine unrestricted_hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight,
Fc(:,:,:) = 0d0
case('LYP')
case('B3LYP')
allocate(FcLDA(nBas,nBas,nspin),FcGGA(nBas,nBas,nspin))
@ -47,13 +47,17 @@ subroutine unrestricted_hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight,
Fc(:,:,:) = FcLDA(:,:,:) + aC*(FcGGA(:,:,:) - FcLDA(:,:,:))
case('BHHLYP')
allocate(FcGGA(nBas,nBas,nspin))
call unrestricted_gga_correlation_potential('LYP ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
case('PBE')
allocate(FcGGA(nBas,nBas,nspin))
call unrestricted_gga_correlation_potential('PBE ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FcGGA)
Fc(:,:,:) = FcGGA(:,:,:)
call unrestricted_gga_correlation_potential('PBE ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
case default

9
src/eDFT/unrestricted_hybrid_exchange_energy.f90
View File

@ -38,7 +38,7 @@ subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1
call unrestricted_fock_exchange_energy(nBas,P,FxHF,Ex)
case ('B3')
case ('B3LYP')
a0 = 0.20d0
aX = 0.72d0
@ -52,6 +52,13 @@ subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1
+ a0*(ExHF - ExLDA) &
+ aX*(ExGGA - ExLDA)
case ('BHHLYP')
call unrestricted_gga_exchange_energy('B88 ',nEns,wEns,nGrid,weight,rho,drho,ExGGA)
call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF)
Ex = 0.5d0*ExHF + 0.5d0*ExGGA
case ('PBE')
call unrestricted_gga_exchange_energy('PBE ',nEns,wEns,nGrid,weight,rho,drho,ExGGA)

11
src/eDFT/unrestricted_hybrid_exchange_potential.f90
View File

@ -43,7 +43,7 @@ subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC
call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF)
Fx(:,:) = FxHF(:,:)
case('B3')
case('B3LYP')
allocate(FxLDA(nBas,nBas),FxGGA(nBas,nBas))
@ -59,6 +59,15 @@ subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC
+ a0*(FxHF(:,:) - FxLDA(:,:)) &
+ aX*(FxGGA(:,:) - FxLDA(:,:))
case('BHHLYP')
allocate(FxGGA(nBas,nBas))
call unrestricted_gga_exchange_potential('B88 ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA)
call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF)
Fx(:,:) = 0.5d0*FxHF(:,:) + 0.5d0*FxGGA(:,:)
case('PBE')
allocate(FxGGA(nBas,nBas))

40
src/eDFT/xc_potential.f90 Normal file
View File

@ -0,0 +1,40 @@
subroutine xc_potential(nBas,c,Fx,Fc,Vxc)
! Compute the exchange-correlation potential in the MO basis
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: Fx(nBas,nBas,nspin)
double precision,intent(in) :: Fc(nBas,nBas,nspin)
! Local variables
integer :: mu,nu
integer :: p
integer :: ispin
! Output variables
double precision,intent(out) :: Vxc(nBas,nspin)
! Compute Vxc
Vxc(:,:) = 0d0
do p=1,nBas
do ispin=1,nspin
do mu=1,nBas
do nu=1,nBas
Vxc(p,ispin) = Vxc(p,ispin) &
+ c(mu,p,ispin)*(Fx(mu,nu,ispin) + Fc(mu,nu,ispin))*c(nu,p,ispin)
end do
end do
end do
end do
end subroutine xc_potential