10
1
mirror of https://github.com/pfloos/quack synced 2024-07-04 02:16:06 +02:00

added choice of exchange coefficient and choice of N-centered in the input file

This commit is contained in:
Clotilde Marut 2020-08-02 13:09:30 +02:00
parent e068aab977
commit a4870ed165
38 changed files with 244 additions and 116 deletions

View File

@ -20,15 +20,19 @@
3
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.00 1.00
# Ncentered ? 0 for NO
0
# Parameters for CC weight-dependent exchange functional
-0.00201219 -0.00371002 -0.00212719
-0.00117715 0.00188738 -0.000414761
0.420431 0.069097 -0.295049
0.135075 -0.00770826 -0.028057
# choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2
2
# occupation numbers of orbitals nO and nO+1
1.00 0.00
1.00 0.00
0.00 0.00
1.00 0.00
1.00 0.00
1.00 1.00
1.00 0.00
0.00 0.00
1.00 0.00
1.00 0.00
1.00 1.00
# GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type
1000 0.00001 T 5 1 1

View File

@ -1,5 +1,5 @@
subroutine GOK_RKS(restart,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,nV,S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum,Cx_choice)
! Perform restricted Kohn-Sham calculation for ensembles
@ -37,7 +37,8 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_
double precision,intent(in) :: ENuc
double precision,intent(inout):: c(nBas,nBas)
double precision,intent(inout),dimension(2,2,3):: occnum
double precision,intent(in):: occnum(2,2,32,2,32,2,3)
integer,intent(in) :: Cx_choice
! Local variables
@ -238,7 +239,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_
! Compute exchange potential
call restricted_exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,Pw(:,:), &
ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fx(:,:),FxHF(:,:))
ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fx(:,:),FxHF(:,:),Cx_choice)
! Compute correlation potential
@ -296,7 +297,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_
! Exchange energy
call restricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, &
Pw(:,:),FxHF(:,:),rhow(:),drhow(:,:),Ex)
Pw(:,:),FxHF(:,:),rhow(:),drhow(:,:),Ex,Cx_choice)
! Correlation energy
@ -345,6 +346,6 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_
call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:), &
nBas,nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc,eps(:),Pw(:,:),rhow(:),drhow(:,:), &
J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,E(:),Om(:),occnum)
J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,E(:),Om(:),occnum,Cx_choice)
end subroutine GOK_RKS

View File

@ -1,5 +1,5 @@
subroutine 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)
nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice)
! Perform unrestricted Kohn-Sham calculation for ensembles
@ -30,7 +30,8 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(in),dimension(2,2,3) :: occnum
double precision,intent(in) :: occnum(2,2,3)
integer,intent(in) :: Cx_choice
! Local variables
@ -239,7 +240,7 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
do ispin=1,nspin
call unrestricted_exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,&
Pw(:,:,ispin),ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin), &
Fx(:,:,ispin),FxHF(:,:,ispin))
Fx(:,:,ispin),FxHF(:,:,ispin),Cx_choice)
end do
! Compute correlation potential
@ -317,7 +318,7 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
do ispin=1,nspin
call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin))
Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin),Cx_choice)
end do
! Correlation energy
@ -369,6 +370,6 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
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,nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum)
AO,dAO,nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum,Cx_choice)
end subroutine GOK_UKS

View File

@ -1,5 +1,6 @@
subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGrid,weight, &
maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c,occnum)
maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,&
c,occnum,Cx_choice)
! Perform restricted Kohn-Sham calculation for ensembles
@ -30,7 +31,8 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(in),dimension(2,2,3) :: occnum
double precision,intent(in) :: occnum(2,2,3)
integer,intent(in) :: Cx_choice
double precision,intent(out) :: c(nBas,nBas)
@ -75,7 +77,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
write(*,*)
call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c,occnum)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c,occnum,Cx_choice)
!------------------------------------------------------------------------
! Equiensemble calculation
@ -96,7 +98,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
write(*,*)
call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c,occnum)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c,occnum,Cx_choice)
!------------------------------------------------------------------------
! Equiensemble calculation
@ -117,7 +119,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
write(*,*)
! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c,occnum)
! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c,occnum,Cx_choice)
!------------------------------------------------------------------------

View File

@ -1,5 +1,6 @@
subroutine 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,nV,S,T,V,Hc,ERI,X,ENuc,c,occnum)
aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,&
nO,nV,S,T,V,Hc,ERI,X,ENuc,c,occnum,Cx_choice)
! Perform restricted Kohn-Sham calculation for ensembles
@ -30,7 +31,8 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(in),dimension(2,2,3) :: occnum
double precision,intent(in) :: occnum(2,2,3)
integer,intent(in) :: Cx_choice
double precision,intent(out) :: c(nBas,nBas)
@ -75,7 +77,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
write(*,*)
call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c,occnum)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c,occnum,Cx_choice)
!------------------------------------------------------------------------
! Equiensemble calculation
@ -96,7 +98,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
write(*,*)
! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c,occnum)
! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c,occnum,Cx_choice)
!------------------------------------------------------------------------
! Equiensemble calculation
@ -117,7 +119,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
write(*,*)
call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c,occnum)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c,occnum,Cx_choice)
!------------------------------------------------------------------------
! MOM excitation energies

View File

@ -1,4 +1,4 @@
subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD)
subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD,Cx_choice)
! Compute the restricted version of the curvature-corrected exchange ensemble derivative
@ -14,6 +14,7 @@ subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -70,15 +71,25 @@ subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
w1 = wEns(2)
w2 = wEns(3)
dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) &
* (1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2))
select case (Cx_choice)
case(1)
dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0)))
dCxdw2 = 0.d0
case(2)
dCxdw1 = 0.d0
dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
case(3)
dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) &
* (1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2))
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)))
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)))
end select
dCxdw1 = CxLDA*dCxdw1
dCxdw2 = CxLDA*dCxdw2
dExdw(:) = 0d0
do iG=1,nGrid

View File

@ -1,4 +1,4 @@
subroutine RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
subroutine RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,Cx_choice)
! Compute the restricted version of the curvature-corrected exchange functional
@ -14,6 +14,7 @@ subroutine RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -67,7 +68,14 @@ subroutine RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
w2 = wEns(3)
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
Cx = CxLDA*Fx1*Fx2
select case (Cx_choice)
case(1)
Cx = CxLDA*Fx1
case(2)
Cx = CxLDA*Fx2
case(3)
Cx = CxLDA*Fx2*Fx1
end select
! Compute GIC-LDA exchange energy

View File

@ -1,4 +1,4 @@
subroutine RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex)
subroutine RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
! Compute the restricted version of the curvature-corrected exchange functional
@ -15,6 +15,7 @@ subroutine RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -73,7 +74,14 @@ subroutine RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
w2 = wEns(3)
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
Cx = CxLDA*Fx1*Fx2
select case (Cx_choice)
case(1)
Cx = CxLDA*Fx1
case(2)
Cx = CxLDA*Fx2
case(3)
Cx = CxLDA*Fx2*Fx1
end select
! Compute LDA exchange matrix in the AO basis

View File

@ -1,4 +1,4 @@
subroutine RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx)
subroutine RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice)
! Compute the restricted version of the curvature-corrected exchange potential
@ -16,6 +16,7 @@ subroutine RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -72,7 +73,15 @@ subroutine RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
w2 = wEns(3)
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
Cx = CxLDA*Fx1*Fx2
select case (Cx_choice)
case(1)
Cx = CxLDA*Fx1
case(2)
Cx = CxLDA*Fx2
case(3)
Cx = CxLDA*Fx2*Fx1
end select
! Compute LDA exchange matrix in the AO basis

View File

@ -1,4 +1,4 @@
subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD)
subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD,Cx_choice)
! Compute the unrestricted version of the curvature-corrected exchange ensemble derivative
@ -14,6 +14,7 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -80,6 +81,21 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
w1 = wEns(2)
w2 = wEns(3)
select case (Cx_choice)
case(1)
dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0)))
dCxdw2 = 0.d0
case(2)
dCxdw1 = 0.d0
dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
case(3)
dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) &
* (1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2))
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)))
end select
! Double weight-dependency
! dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) &
@ -93,8 +109,8 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
! dCxdw2 = 0.d0
! right single-weight-dependency
dCxdw1 = 0.d0
dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
! dCxdw1 = 0.d0
! dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))

View File

@ -1,4 +1,4 @@
subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,Cx_choice)
! Compute the unrestricted version of the curvature-corrected exchange functional
@ -14,6 +14,7 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -77,6 +78,14 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
w2 = wEns(3)
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
end select
! for two-weights ensemble
! Cx = alpha*Fx2*Fx1
@ -84,7 +93,7 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
! Cx = alpha*Fx1
! for right ensemble
Cx = alpha*Fx2
! Cx = alpha*Fx2
! Compute GIC-LDA exchange energy

View File

@ -1,4 +1,4 @@
subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex)
subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
! Compute the unrestricted version of the curvature-corrected exchange functional
@ -15,7 +15,8 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Local variables
integer :: iG
@ -76,6 +77,15 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
w2 = wEns(3)
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
end select
! for two-weight ensembles
! Cx = alpha*Fx1*Fx2
@ -83,7 +93,7 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
! Cx = alpha*Fx1
! for right ensembles
Cx = alpha*Fx2
! Cx = alpha*Fx2
! Compute LDA exchange matrix in the AO basis

View File

@ -1,4 +1,4 @@
subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx)
subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice)
! Compute the unrestricted version of the curvature-corrected exchange potential
@ -16,6 +16,7 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -79,6 +80,15 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
w2 = wEns(3)
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
end select
! for two-weight ensembles
! Cx = alpha*Fx2*Fx1
@ -86,7 +96,7 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
! Cx = alpha*Fx1
! for right ensemble
Cx = alpha*Fx2
! Cx = alpha*Fx2
! Compute LDA exchange matrix in the AO basis

View File

@ -50,7 +50,7 @@ program eDFT
double precision :: start_KS,end_KS,t_KS
double precision :: start_int,end_int,t_int
integer :: nEns
integer :: nEns,doNcentered
double precision,allocatable :: wEns(:)
integer :: maxSCF,max_diis
@ -59,7 +59,9 @@ program eDFT
integer :: guess_type
integer :: ortho_type
double precision,dimension(2,2,3) :: occnum
! double precision,allocatable,dimension(:,:,:) :: occnum
double precision,dimension(2,2,3) :: occnum
integer :: Cx_choice
! Hello World
@ -113,7 +115,7 @@ program eDFT
allocate(wEns(maxEns))
call read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, &
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,ncent,occnum)
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,doNcentered,ncent,occnum,Cx_choice)
!------------------------------------------------------------------------
! Read one- and two-electron integrals
@ -175,7 +177,7 @@ program eDFT
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)
S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum,Cx_choice)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -193,7 +195,7 @@ program eDFT
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)
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -211,7 +213,7 @@ program eDFT
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)
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -228,7 +230,7 @@ program eDFT
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)
nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum,Cx_choice)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -245,7 +247,7 @@ program eDFT
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, &
nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum)
nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice)
call cpu_time(end_KS)
t_KS = end_KS - start_KS

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, &
nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum)
nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice)
! Perform unrestricted Kohn-Sham calculation for ensembles
@ -30,8 +30,8 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(in),dimension(2,2,3) :: occnum
double precision,intent(in) :: occnum(2,2,3)
integer,intent(in) :: Cx_choice
! Local variables
@ -242,7 +242,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
do ispin=1,nspin
call unrestricted_exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, &
Pw(:,:,ispin),ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin), &
Fx(:,:,ispin),FxHF(:,:,ispin))
Fx(:,:,ispin),FxHF(:,:,ispin),Cx_choice)
end do
! Compute correlation potential
@ -320,7 +320,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
do ispin=1,nspin
call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin))
Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin),Cx_choice)
end do
! Correlation energy
@ -371,6 +371,6 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
!------------------------------------------------------------------------
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,nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum)
AO,dAO,nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum,Cx_choice)
end subroutine eDFT_UKS

View File

@ -1,5 +1,5 @@
subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, &
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,ncent,occnum)
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,doNcentered,ncent,occnum,Cx_choice)
! Read DFT options
@ -17,11 +17,11 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
integer,intent(out) :: x_rung,c_rung
character(len=12),intent(out) :: x_DFA, c_DFA
integer,intent(out) :: SGn
integer,intent(out) :: nEns
integer,intent(out) :: nEns, doNcentered
double precision,intent(out) :: wEns(maxEns)
double precision,intent(out) :: aCC_w1(3)
double precision,intent(out) :: aCC_w2(3)
double precision,intent(out),dimension(2,2,3) :: occnum
double precision,intent(inout) :: occnum(2,2,3)
integer,intent(out) :: maxSCF
double precision,intent(out) :: thresh
@ -30,6 +30,7 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
integer,intent(out) :: guess_type
integer,intent(out) :: ortho_type
double precision,intent(in) :: ncent
integer,intent(out) :: Cx_choice
! Local variables
@ -105,8 +106,15 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
! Read ensemble weights for real physical (fractional number of electrons) ensemble (w1,w2)
read(1,*)
read(1,*) (wEns(I),I=2,nEns)
wEns(1) = 1d0 - wEns(2) - wEns(3)
read(1,*)
read(1,*) doNcentered
if (doNcentered==0) then
wEns(1) = 1d0 - wEns(2) - wEns(3)
else
wEns(2) = (ncent/(ncent-1.d0))*wEns(2)
wEns(3) = (ncent/(ncent+1.d0))*wEns(3)
wEns(1) = 1d0 - ((ncent-1.d0)/ncent)*wEns(2) - ((ncent+1.d0)/ncent)*wEns(3) ! for N-centered
end if
write(*,*)'----------------------------------------------------------'
write(*,*)' Ensemble weights '
write(*,*)'----------------------------------------------------------'
@ -117,6 +125,9 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
read(1,*)
read(1,*) (aCC_w1(I),I=1,3)
read(1,*) (aCC_w2(I),I=1,3)
! Read choice of exchange coefficient
read(1,*)
read(1,*) Cx_choice
write(*,*)'----------------------------------------------------------'
write(*,*)' parameters for w1-dependant exchange functional coefficient '
@ -130,7 +141,7 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
call matout(3,1,aCC_w2)
write(*,*)
! allocate(occnum(2,2,nEns))
!allocate(occnum(nspin,2,nEns))
! Read occupation numbers for orbitals nO and nO+1
read(1,*)
do J=1,3

View File

@ -12,7 +12,7 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
integer,intent(in) :: nEns
integer,intent(in) :: nO
double precision,intent(in) :: eps(nBas)
double precision,intent(in),dimension(2,2,3) :: occnum
double precision,intent(in) :: occnum(2,2,3)
! Local variables

View File

@ -12,7 +12,7 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P,occnum)
integer,intent(in) :: nEns
integer,intent(in) :: nO
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in),dimension(2,2,3) :: occnum
double precision,intent(in) :: occnum(2,2,3)
! Local variables

View File

@ -1,4 +1,4 @@
subroutine restricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,ExDD)
subroutine restricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,ExDD,Cx_choice)
! Compute the exchange part of the derivative discontinuity
@ -17,6 +17,7 @@ subroutine restricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -37,7 +38,8 @@ subroutine restricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w
case(1)
call restricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:))
call restricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),&
rhow(:),ExDD(:),Cx_choice)
! GGA functionals

View File

@ -1,4 +1,4 @@
subroutine restricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF,rho,drho,Ex)
subroutine restricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF,rho,drho,Ex,Cx_choice)
! Compute the exchange energy
@ -21,6 +21,7 @@ subroutine restricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC
double precision,intent(in) :: FxHF(nBas,nBas)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -43,7 +44,7 @@ subroutine restricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC
case(1)
call restricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,ExLDA)
call restricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,ExLDA,Cx_choice)
Ex = ExLDA
@ -63,7 +64,7 @@ subroutine restricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC
aX = 0.72d0
aC = 0.81d0
call restricted_lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,ExLDA)
call restricted_lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,ExLDA,Cx_choice)
call restricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA)
call restricted_fock_exchange_energy(nBas,P,FxHF,ExHF)

View File

@ -1,5 +1,5 @@
subroutine restricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
ERI,Pw,P,rhow,drhow,rho,drho,Ex)
ERI,Pw,P,rhow,drhow,rho,drho,Ex,Cx_choice)
! Compute the exchange individual energy
@ -25,6 +25,7 @@ subroutine restricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns
double precision,intent(in) :: drhow(ncart,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -48,7 +49,7 @@ subroutine restricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns
case(1)
call restricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,ExLDA)
call restricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,ExLDA,Cx_choice)
Ex = ExLDA

View File

@ -1,5 +1,5 @@
subroutine restricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P, &
ERI,AO,dAO,rho,drho,Fx,FxHF)
ERI,AO,dAO,rho,drho,Fx,FxHF,Cx_choice)
! Compute the exchange potential
@ -24,6 +24,7 @@ subroutine restricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -48,7 +49,7 @@ subroutine restricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,
case(1)
call restricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx)
call restricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice)
! GGA functionals
@ -65,7 +66,7 @@ subroutine restricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,
cX = 0.20d0
aX = 0.72d0
call restricted_lda_exchange_potential(DFA,nGrid,weight,nBas,AO,rho,FxLDA)
call restricted_lda_exchange_potential(DFA,nGrid,weight,nBas,AO,rho,FxLDA,Cx_choice)
call restricted_gga_exchange_potential(DFA,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA)
call restricted_fock_exchange_potential(nBas,P,ERI,FxHF)

View File

@ -1,5 +1,5 @@
subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,E,Om,occnum)
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,E,Om,occnum,Cx_choice)
! Compute individual energies as well as excitation energies
@ -38,7 +38,8 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n
double precision,intent(in) :: J(nBas,nBas)
double precision :: Ew
double precision,intent(in),dimension(2,2,3) :: occnum
double precision,intent(in) :: occnum(2,2,3)
integer,intent(in) :: Cx_choice
! Local variables
@ -94,7 +95,8 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n
do iEns=1,nEns
call restricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, &
ERI(:,:,:,:),Pw(:,:),P(:,:,iEns),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns),Ex(iEns))
ERI(:,:,:,:),Pw(:,:),P(:,:,iEns),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns)&
,Ex(iEns),Cx_choice)
end do
!------------------------------------------------------------------------
@ -117,7 +119,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n
!------------------------------------------------------------------------
call restricted_exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:), &
rhow(:),drhow(:,:),ExDD(:))
rhow(:),drhow(:,:),ExDD(:),Cx_choice)
call restricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),EcDD(:))

View File

@ -1,4 +1,4 @@
subroutine restricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD)
subroutine restricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD,Cx_choice)
! Compute the exchange LDA part of the derivative discontinuity
@ -16,6 +16,7 @@ subroutine restricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -38,7 +39,7 @@ subroutine restricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1
case ('CC')
call RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:))
call RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:),Cx_choice)
case default

View File

@ -1,4 +1,4 @@
subroutine restricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
subroutine restricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,Cx_choice)
! Select LDA exchange functional
@ -16,6 +16,7 @@ subroutine restricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Output variables
@ -35,7 +36,7 @@ subroutine restricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_
case ('CC')
call RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
call RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,Cx_choice)
case default

View File

@ -1,4 +1,4 @@
subroutine restricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex)
subroutine restricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
! Compute LDA exchange energy for individual states
@ -17,6 +17,7 @@ subroutine restricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Output variables
@ -36,7 +37,7 @@ subroutine restricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,
case ('CC')
call RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex)
call RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
case default

View File

@ -1,4 +1,4 @@
subroutine restricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx)
subroutine restricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice)
! Select LDA correlation potential
@ -19,6 +19,7 @@ subroutine restricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,a
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Output variables
@ -38,7 +39,7 @@ subroutine restricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,a
case ('CC')
call RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx)
call RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice)
case default

View File

@ -12,7 +12,7 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
integer,intent(in) :: nEns
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: eps(nBas,nspin)
double precision,intent(in),dimension(2,2,3) :: occnum
double precision,intent(in) :: occnum(2,2,3)
! Local variables

View File

@ -12,7 +12,7 @@ subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P,occnum)
integer,intent(in) :: nEns
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in),dimension(2,2,3) :: occnum
double precision,intent(in) :: occnum(2,2,3)
! Local variables

View File

@ -1,4 +1,4 @@
subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,ExDD)
subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,ExDD,Cx_choice)
! Compute the exchange part of the derivative discontinuity
@ -17,6 +17,7 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -37,7 +38,8 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC
case(1)
call unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:))
call unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),&
rhow(:),ExDD(:),Cx_choice)
! GGA functionals

View File

@ -1,4 +1,5 @@
subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF,rho,drho,Ex)
subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF, &
rho,drho,Ex,Cx_choice)
! Compute the exchange energy
@ -21,6 +22,7 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,a
double precision,intent(in) :: FxHF(nBas,nBas)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -43,7 +45,8 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,a
case(1)
call unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,ExLDA)
call unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,&
rho,ExLDA,Cx_choice)
Ex = ExLDA
@ -63,7 +66,7 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,a
aX = 0.72d0
aC = 0.81d0
call unrestricted_lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,ExLDA)
call unrestricted_lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,ExLDA,Cx_choice)
call unrestricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA)
call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF)

View File

@ -1,5 +1,5 @@
subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
ERI,Pw,P,rhow,drhow,rho,drho,Ex)
ERI,Pw,P,rhow,drhow,rho,drho,Ex,Cx_choice)
! Compute the exchange individual energy
@ -25,6 +25,7 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE
double precision,intent(in) :: drhow(ncart,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -48,7 +49,8 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE
case(1)
call unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,ExLDA)
call unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,&
rhow,rho,ExLDA,Cx_choice)
Ex = ExLDA

View File

@ -1,5 +1,5 @@
subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P, &
ERI,AO,dAO,rho,drho,Fx,FxHF)
ERI,AO,dAO,rho,drho,Fx,FxHF,Cx_choice)
! Compute the exchange potential
@ -24,6 +24,7 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -48,7 +49,7 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w
case(1)
call unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx)
call unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice)
! GGA functionals
@ -65,7 +66,7 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w
cX = 0.20d0
aX = 0.72d0
call unrestricted_lda_exchange_potential(DFA,nGrid,weight,nBas,AO,rho,FxLDA)
call unrestricted_lda_exchange_potential(DFA,nGrid,weight,nBas,AO,rho,FxLDA,Cx_choice)
call unrestricted_gga_exchange_potential(DFA,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA)
call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF)

View File

@ -1,5 +1,5 @@
subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,dAO, &
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum)
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum,Cx_choice)
! Compute unrestricted individual energies as well as excitation energies
@ -41,7 +41,8 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
double precision,intent(in) :: FxHF(nBas,nBas,nspin)
double precision,intent(in) :: Fc(nBas,nBas,nspin)
double precision :: Ew
double precision,intent(in),dimension(2,2,3) :: occnum
double precision,intent(in) :: occnum(2,2,3)
integer,intent(in) :: Cx_choice
! Local variables
@ -139,7 +140,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
do ispin=1,nspin
call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,ERI, &
Pw(:,:,ispin),P(:,:,ispin,iEns),rhow(:,ispin),drhow(:,:,ispin), &
rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns))
rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns),Cx_choice)
end do
end do
@ -193,7 +194,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
do ispin=1,nspin
call unrestricted_exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, &
rhow(:,ispin),drhow(:,:,ispin),ExDD(ispin,:))
rhow(:,ispin),drhow(:,:,ispin),ExDD(ispin,:),Cx_choice)
end do
call unrestricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD)

View File

@ -1,4 +1,4 @@
subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD)
subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD,Cx_choice)
! Compute the exchange LDA part of the derivative discontinuity
@ -16,6 +16,7 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
integer,intent(in) :: Cx_choice
! Local variables
@ -34,7 +35,7 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_
case ('CC')
call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:))
call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:),Cx_choice)
case default

View File

@ -1,4 +1,4 @@
subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,Cx_choice)
! Select LDA exchange functional
@ -16,6 +16,7 @@ subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aC
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Output variables
@ -31,7 +32,7 @@ subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aC
case ('CC')
call UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
call UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,Cx_choice)
case default

View File

@ -1,4 +1,4 @@
subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex)
subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
! Compute LDA exchange energy for individual states
@ -17,6 +17,7 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Output variables
@ -32,7 +33,7 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn
case ('CC')
call UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex)
call UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
case default

View File

@ -1,4 +1,4 @@
subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx)
subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice)
! Select LDA correlation potential
@ -19,6 +19,7 @@ subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Output variables
@ -34,7 +35,7 @@ subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1
case ('CC')
call UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx)
call UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice)
case default