From a4870ed1650b30572c5e9a19e718cb0ad1b67874 Mon Sep 17 00:00:00 2001 From: Clotilde Marut Date: Sun, 2 Aug 2020 13:09:30 +0200 Subject: [PATCH] added choice of exchange coefficient and choice of N-centered in the input file --- input/dft | 18 +++++++++------ src/eDFT/GOK_RKS.f90 | 11 +++++---- src/eDFT/GOK_UKS.f90 | 11 +++++---- src/eDFT/LIM_RKS.f90 | 12 ++++++---- src/eDFT/MOM_RKS.f90 | 12 ++++++---- ..._lda_exchange_derivative_discontinuity.f90 | 21 +++++++++++++---- src/eDFT/RCC_lda_exchange_energy.f90 | 12 ++++++++-- .../RCC_lda_exchange_individual_energy.f90 | 12 ++++++++-- src/eDFT/RCC_lda_exchange_potential.f90 | 13 +++++++++-- ..._lda_exchange_derivative_discontinuity.f90 | 22 +++++++++++++++--- src/eDFT/UCC_lda_exchange_energy.f90 | 13 +++++++++-- .../UCC_lda_exchange_individual_energy.f90 | 16 ++++++++++--- src/eDFT/UCC_lda_exchange_potential.f90 | 14 +++++++++-- src/eDFT/eDFT.f90 | 18 ++++++++------- src/eDFT/eDFT_UKS.f90 | 12 +++++----- src/eDFT/read_options.f90 | 23 ++++++++++++++----- src/eDFT/restricted_auxiliary_energy.f90 | 2 +- src/eDFT/restricted_density_matrix.f90 | 2 +- ...cted_exchange_derivative_discontinuity.f90 | 6 +++-- src/eDFT/restricted_exchange_energy.f90 | 7 +++--- .../restricted_exchange_individual_energy.f90 | 5 ++-- src/eDFT/restricted_exchange_potential.f90 | 7 +++--- src/eDFT/restricted_individual_energy.f90 | 10 ++++---- ..._lda_exchange_derivative_discontinuity.f90 | 5 ++-- src/eDFT/restricted_lda_exchange_energy.f90 | 5 ++-- ...tricted_lda_exchange_individual_energy.f90 | 5 ++-- .../restricted_lda_exchange_potential.f90 | 5 ++-- src/eDFT/unrestricted_auxiliary_energy.f90 | 2 +- src/eDFT/unrestricted_density_matrix.f90 | 2 +- ...cted_exchange_derivative_discontinuity.f90 | 6 +++-- src/eDFT/unrestricted_exchange_energy.f90 | 9 +++++--- ...nrestricted_exchange_individual_energy.f90 | 6 +++-- src/eDFT/unrestricted_exchange_potential.f90 | 7 +++--- src/eDFT/unrestricted_individual_energy.f90 | 9 ++++---- ..._lda_exchange_derivative_discontinuity.f90 | 5 ++-- src/eDFT/unrestricted_lda_exchange_energy.f90 | 5 ++-- ...tricted_lda_exchange_individual_energy.f90 | 5 ++-- .../unrestricted_lda_exchange_potential.f90 | 5 ++-- 38 files changed, 244 insertions(+), 116 deletions(-) diff --git a/input/dft b/input/dft index 985b124..7ea834b 100644 --- a/input/dft +++ b/input/dft @@ -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 diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index d1c32f7..17b3758 100644 --- a/src/eDFT/GOK_RKS.f90 +++ b/src/eDFT/GOK_RKS.f90 @@ -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 diff --git a/src/eDFT/GOK_UKS.f90 b/src/eDFT/GOK_UKS.f90 index 1fd9ee4..74d91da 100644 --- a/src/eDFT/GOK_UKS.f90 +++ b/src/eDFT/GOK_UKS.f90 @@ -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 diff --git a/src/eDFT/LIM_RKS.f90 b/src/eDFT/LIM_RKS.f90 index 85b8f70..cbaed5b 100644 --- a/src/eDFT/LIM_RKS.f90 +++ b/src/eDFT/LIM_RKS.f90 @@ -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) !------------------------------------------------------------------------ diff --git a/src/eDFT/MOM_RKS.f90 b/src/eDFT/MOM_RKS.f90 index cd7e00e..1f6ec10 100644 --- a/src/eDFT/MOM_RKS.f90 +++ b/src/eDFT/MOM_RKS.f90 @@ -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 diff --git a/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 index 56370c6..235407a 100644 --- a/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 @@ -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 diff --git a/src/eDFT/RCC_lda_exchange_energy.f90 b/src/eDFT/RCC_lda_exchange_energy.f90 index 7a6d1a5..8fa5c6c 100644 --- a/src/eDFT/RCC_lda_exchange_energy.f90 +++ b/src/eDFT/RCC_lda_exchange_energy.f90 @@ -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 diff --git a/src/eDFT/RCC_lda_exchange_individual_energy.f90 b/src/eDFT/RCC_lda_exchange_individual_energy.f90 index 66a272e..02f23f8 100644 --- a/src/eDFT/RCC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RCC_lda_exchange_individual_energy.f90 @@ -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 diff --git a/src/eDFT/RCC_lda_exchange_potential.f90 b/src/eDFT/RCC_lda_exchange_potential.f90 index a75bdfe..d33017b 100644 --- a/src/eDFT/RCC_lda_exchange_potential.f90 +++ b/src/eDFT/RCC_lda_exchange_potential.f90 @@ -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 diff --git a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 index caf6122..1246d3b 100644 --- a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 @@ -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))) diff --git a/src/eDFT/UCC_lda_exchange_energy.f90 b/src/eDFT/UCC_lda_exchange_energy.f90 index 3618298..b3fcbaf 100644 --- a/src/eDFT/UCC_lda_exchange_energy.f90 +++ b/src/eDFT/UCC_lda_exchange_energy.f90 @@ -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 diff --git a/src/eDFT/UCC_lda_exchange_individual_energy.f90 b/src/eDFT/UCC_lda_exchange_individual_energy.f90 index 25fe2cd..72867b3 100644 --- a/src/eDFT/UCC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/UCC_lda_exchange_individual_energy.f90 @@ -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 diff --git a/src/eDFT/UCC_lda_exchange_potential.f90 b/src/eDFT/UCC_lda_exchange_potential.f90 index 2dce16a..dc14909 100644 --- a/src/eDFT/UCC_lda_exchange_potential.f90 +++ b/src/eDFT/UCC_lda_exchange_potential.f90 @@ -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 diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 857c733..bc2fbfc 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -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 diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index dd66abd..8ba1589 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -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 diff --git a/src/eDFT/read_options.f90 b/src/eDFT/read_options.f90 index e92eb7f..884ac4b 100644 --- a/src/eDFT/read_options.f90 +++ b/src/eDFT/read_options.f90 @@ -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 diff --git a/src/eDFT/restricted_auxiliary_energy.f90 b/src/eDFT/restricted_auxiliary_energy.f90 index 8c294e0..7bbe751 100644 --- a/src/eDFT/restricted_auxiliary_energy.f90 +++ b/src/eDFT/restricted_auxiliary_energy.f90 @@ -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 diff --git a/src/eDFT/restricted_density_matrix.f90 b/src/eDFT/restricted_density_matrix.f90 index 661113e..d96f7e6 100644 --- a/src/eDFT/restricted_density_matrix.f90 +++ b/src/eDFT/restricted_density_matrix.f90 @@ -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 diff --git a/src/eDFT/restricted_exchange_derivative_discontinuity.f90 b/src/eDFT/restricted_exchange_derivative_discontinuity.f90 index 5b7441b..1d9561a 100644 --- a/src/eDFT/restricted_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/restricted_exchange_derivative_discontinuity.f90 @@ -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 diff --git a/src/eDFT/restricted_exchange_energy.f90 b/src/eDFT/restricted_exchange_energy.f90 index c9e1d67..fc43d4e 100644 --- a/src/eDFT/restricted_exchange_energy.f90 +++ b/src/eDFT/restricted_exchange_energy.f90 @@ -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) diff --git a/src/eDFT/restricted_exchange_individual_energy.f90 b/src/eDFT/restricted_exchange_individual_energy.f90 index a4aa93f..12c787a 100644 --- a/src/eDFT/restricted_exchange_individual_energy.f90 +++ b/src/eDFT/restricted_exchange_individual_energy.f90 @@ -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 diff --git a/src/eDFT/restricted_exchange_potential.f90 b/src/eDFT/restricted_exchange_potential.f90 index 81fcea0..ddcbb2b 100644 --- a/src/eDFT/restricted_exchange_potential.f90 +++ b/src/eDFT/restricted_exchange_potential.f90 @@ -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) diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index 466e300..1d4cd92 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -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(:)) diff --git a/src/eDFT/restricted_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/restricted_lda_exchange_derivative_discontinuity.f90 index 10eebb9..8223222 100644 --- a/src/eDFT/restricted_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/restricted_lda_exchange_derivative_discontinuity.f90 @@ -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 diff --git a/src/eDFT/restricted_lda_exchange_energy.f90 b/src/eDFT/restricted_lda_exchange_energy.f90 index 519c0a8..ae6d75c 100644 --- a/src/eDFT/restricted_lda_exchange_energy.f90 +++ b/src/eDFT/restricted_lda_exchange_energy.f90 @@ -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 diff --git a/src/eDFT/restricted_lda_exchange_individual_energy.f90 b/src/eDFT/restricted_lda_exchange_individual_energy.f90 index fcdde39..c72c181 100644 --- a/src/eDFT/restricted_lda_exchange_individual_energy.f90 +++ b/src/eDFT/restricted_lda_exchange_individual_energy.f90 @@ -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 diff --git a/src/eDFT/restricted_lda_exchange_potential.f90 b/src/eDFT/restricted_lda_exchange_potential.f90 index 9f1b5bf..e2f740b 100644 --- a/src/eDFT/restricted_lda_exchange_potential.f90 +++ b/src/eDFT/restricted_lda_exchange_potential.f90 @@ -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 diff --git a/src/eDFT/unrestricted_auxiliary_energy.f90 b/src/eDFT/unrestricted_auxiliary_energy.f90 index be14e03..cb6e680 100644 --- a/src/eDFT/unrestricted_auxiliary_energy.f90 +++ b/src/eDFT/unrestricted_auxiliary_energy.f90 @@ -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 diff --git a/src/eDFT/unrestricted_density_matrix.f90 b/src/eDFT/unrestricted_density_matrix.f90 index 36de196..41b73d1 100644 --- a/src/eDFT/unrestricted_density_matrix.f90 +++ b/src/eDFT/unrestricted_density_matrix.f90 @@ -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 diff --git a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 index 04aae3a..30143bc 100644 --- a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 @@ -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 diff --git a/src/eDFT/unrestricted_exchange_energy.f90 b/src/eDFT/unrestricted_exchange_energy.f90 index f31c2aa..91c8e52 100644 --- a/src/eDFT/unrestricted_exchange_energy.f90 +++ b/src/eDFT/unrestricted_exchange_energy.f90 @@ -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) diff --git a/src/eDFT/unrestricted_exchange_individual_energy.f90 b/src/eDFT/unrestricted_exchange_individual_energy.f90 index 3a0844e..b0d5a66 100644 --- a/src/eDFT/unrestricted_exchange_individual_energy.f90 +++ b/src/eDFT/unrestricted_exchange_individual_energy.f90 @@ -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 diff --git a/src/eDFT/unrestricted_exchange_potential.f90 b/src/eDFT/unrestricted_exchange_potential.f90 index 573a70e..487e5b4 100644 --- a/src/eDFT/unrestricted_exchange_potential.f90 +++ b/src/eDFT/unrestricted_exchange_potential.f90 @@ -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) diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 index 1e6d58d..1b4dd22 100644 --- a/src/eDFT/unrestricted_individual_energy.f90 +++ b/src/eDFT/unrestricted_individual_energy.f90 @@ -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) diff --git a/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 index e84be1a..4750bd1 100644 --- a/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 @@ -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 diff --git a/src/eDFT/unrestricted_lda_exchange_energy.f90 b/src/eDFT/unrestricted_lda_exchange_energy.f90 index 76829b4..f4d6fba 100644 --- a/src/eDFT/unrestricted_lda_exchange_energy.f90 +++ b/src/eDFT/unrestricted_lda_exchange_energy.f90 @@ -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 diff --git a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 index 56b6a4a..9d64dfe 100644 --- a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 +++ b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 @@ -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 diff --git a/src/eDFT/unrestricted_lda_exchange_potential.f90 b/src/eDFT/unrestricted_lda_exchange_potential.f90 index 9f47532..dd26ffb 100644 --- a/src/eDFT/unrestricted_lda_exchange_potential.f90 +++ b/src/eDFT/unrestricted_lda_exchange_potential.f90 @@ -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