From 2eddd491bf2abe273fe2f60ad888e204946e3c38 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 1 Apr 2020 22:59:52 +0200 Subject: [PATCH] individual energies fixed at fucking last --- include/parameters.h | 4 +- input/basis | 54 ++++++++------- input/dft | 4 +- input/molecule | 5 +- input/molecule.xyz | 5 +- input/weight | 54 ++++++++------- src/eDFT/GOK_RKS.f90 | 22 ++++--- src/eDFT/LIM_RKS.f90 | 9 +-- src/eDFT/RMFL20_lda_correlation_energy.f90 | 6 +- ...FL20_lda_correlation_individual_energy.f90 | 5 +- src/eDFT/RMFL20_lda_correlation_potential.f90 | 4 +- src/eDFT/RMFL20_lda_exchange_energy.f90 | 4 +- .../RMFL20_lda_exchange_individual_energy.f90 | 6 +- src/eDFT/RMFL20_lda_exchange_potential.f90 | 4 +- .../RS51_lda_exchange_individual_energy.f90 | 3 +- ...VWN5_lda_correlation_individual_energy.f90 | 12 ++-- src/eDFT/RVWN5_lda_correlation_potential.f90 | 10 +-- src/eDFT/eDFT.f90 | 15 +++-- src/eDFT/exchange_energy.f90 | 5 +- src/eDFT/exchange_individual_energy.f90 | 10 +-- src/eDFT/exchange_potential.f90 | 5 +- src/eDFT/fock_exchange_individual_energy.f90 | 3 +- src/eDFT/hartree_individual_energy.f90 | 65 +++++++++++++++++++ src/eDFT/lda_exchange_energy.f90 | 5 +- src/eDFT/lda_exchange_individual_energy.f90 | 5 +- src/eDFT/lda_exchange_potential.f90 | 5 +- src/eDFT/restricted_correlation_energy.f90 | 5 +- ...stricted_correlation_individual_energy.f90 | 5 +- src/eDFT/restricted_correlation_potential.f90 | 5 +- src/eDFT/restricted_individual_energy.f90 | 21 ++---- .../restricted_lda_correlation_energy.f90 | 5 +- ...cted_lda_correlation_individual_energy.f90 | 5 +- .../restricted_lda_correlation_potential.f90 | 5 +- 33 files changed, 236 insertions(+), 144 deletions(-) create mode 100644 src/eDFT/hartree_individual_energy.f90 diff --git a/include/parameters.h b/include/parameters.h index 850485b..ee9572e 100644 --- a/include/parameters.h +++ b/include/parameters.h @@ -18,6 +18,6 @@ double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0) double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) -! double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) - double precision,parameter :: Cx1 = - 0.904*(176d0/105d0)*(1d0/pi)**(1d0/3d0) +! double precision,parameter :: Cx1 = - 0.913d0*(4d0/3d0)*(1d0/pi)**(1d0/3d0) + double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) diff --git a/input/basis b/input/basis index a1f7a8d..6d7be41 100644 --- a/input/basis +++ b/input/basis @@ -1,27 +1,35 @@ -1 5 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 +1 9 +S 8 + 1 2940.0000000 0.0006800 + 2 441.2000000 0.0052360 + 3 100.5000000 0.0266060 + 4 28.4300000 0.0999930 + 5 9.1690000 0.2697020 + 6 3.1960000 0.4514690 + 7 1.1590000 0.2950740 + 8 0.1811000 0.0125870 +S 8 + 1 2940.0000000 -0.0001230 + 2 441.2000000 -0.0009660 + 3 100.5000000 -0.0048310 + 4 28.4300000 -0.0193140 + 5 9.1690000 -0.0532800 + 6 3.1960000 -0.1207230 + 7 1.1590000 -0.1334350 + 8 0.1811000 0.5307670 S 1 - 1 0.1220000 1.0000000 + 1 0.0589000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.0187700 1.0000000 +P 3 + 1 3.6190000 0.0291110 + 2 0.7110000 0.1693650 + 3 0.1951000 0.5134580 P 1 - 1 0.7270000 1.0000000 + 1 0.0601800 1.0000000 P 1 - 1 0.1410000 1.0000000 -2 5 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 -S 1 - 1 0.1220000 1.0000000 -S 1 - 1 0.0297400 1.0000000 -P 1 - 1 0.7270000 1.0000000 -P 1 - 1 0.1410000 1.0000000 - + 1 0.0085000 1.0000000 +D 1 + 1 0.2380000 1.0000000 +D 1 + 1 0.0740000 1.0000000 diff --git a/input/dft b/input/dft index 49205b9..282688b 100644 --- a/input/dft +++ b/input/dft @@ -1,5 +1,5 @@ # Restricted or unrestricted KS calculation - LIM-RKS + GOK-RKS # exchange rung: # Hartree = 0 # LDA = 1: RS51,RMFL20 @@ -19,6 +19,6 @@ # Number of states in ensemble (nEns) 2 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.25 + 0.5 # GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type 32 0.00001 T 5 1 1 diff --git a/input/molecule b/input/molecule index 8076140..6a6f6d1 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,4 @@ # nAt nEla nElb nCore nRyd - 2 1 1 0 0 + 1 2 2 0 0 # Znuc x y z - H 0.0 0.0 0.0 - H 0.0 0.0 1.4 + Be 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index 6edc99d..8023e37 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,3 @@ - 2 + 1 - H 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 0.7408481486 + Be 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/weight b/input/weight index a1f7a8d..6d7be41 100644 --- a/input/weight +++ b/input/weight @@ -1,27 +1,35 @@ -1 5 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 +1 9 +S 8 + 1 2940.0000000 0.0006800 + 2 441.2000000 0.0052360 + 3 100.5000000 0.0266060 + 4 28.4300000 0.0999930 + 5 9.1690000 0.2697020 + 6 3.1960000 0.4514690 + 7 1.1590000 0.2950740 + 8 0.1811000 0.0125870 +S 8 + 1 2940.0000000 -0.0001230 + 2 441.2000000 -0.0009660 + 3 100.5000000 -0.0048310 + 4 28.4300000 -0.0193140 + 5 9.1690000 -0.0532800 + 6 3.1960000 -0.1207230 + 7 1.1590000 -0.1334350 + 8 0.1811000 0.5307670 S 1 - 1 0.1220000 1.0000000 + 1 0.0589000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.0187700 1.0000000 +P 3 + 1 3.6190000 0.0291110 + 2 0.7110000 0.1693650 + 3 0.1951000 0.5134580 P 1 - 1 0.7270000 1.0000000 + 1 0.0601800 1.0000000 P 1 - 1 0.1410000 1.0000000 -2 5 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 -S 1 - 1 0.1220000 1.0000000 -S 1 - 1 0.0297400 1.0000000 -P 1 - 1 0.7270000 1.0000000 -P 1 - 1 0.1410000 1.0000000 - + 1 0.0085000 1.0000000 +D 1 + 1 0.2380000 1.0000000 +D 1 + 1 0.0740000 1.0000000 diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index bfd8db2..5bd67f9 100644 --- a/src/eDFT/GOK_RKS.f90 +++ b/src/eDFT/GOK_RKS.f90 @@ -1,4 +1,4 @@ -subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh, & +subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,maxSCF,thresh, & max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,EwGIC,c) ! Perform restricted Kohn-Sham calculation for ensembles @@ -11,11 +11,14 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS logical,intent(in) :: restart integer,intent(in) :: x_rung,c_rung character(len=12),intent(in) :: x_DFA,c_DFA + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) - integer,intent(in) :: maxSCF,max_diis,guess_type + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + integer,intent(in) :: guess_type double precision,intent(in) :: thresh integer,intent(in) :: nBas double precision,intent(in) :: AO(nBas,nGrid) @@ -232,12 +235,12 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS ! Compute exchange potential - call exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,Pw(:,:),ERI(:,:,:,:), & + call exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas,Pw(:,:),ERI(:,:,:,:), & AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fx(:,:),FxHF(:,:)) ! Compute correlation potential - call restricted_correlation_potential(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:), & + call restricted_correlation_potential(c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), & nBas,AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fc(:,:)) ! Build Fock operator @@ -290,12 +293,12 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS ! Exchange energy - call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas, & + call exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas, & Pw(:,:),FxHF(:,:),rhow(:),drhow(:,:),Ex) ! Correlation energy - call restricted_correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),Ec) + call restricted_correlation_energy(c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),Ec) ! Total energy @@ -338,9 +341,8 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS ! Compute individual energies from ensemble energy !------------------------------------------------------------------------ - call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:), & - nBas,nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc, & - eps(:),Pw(:,:),rhow(:),drhow(:,:),J(:,:),P(:,:,:), & - rho(:,:),drho(:,:,:),Ew,EwGIC,E(:),Om(:)) + call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), & + nBas,nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc,eps(:),Pw(:,:),rhow(:),drhow(:,:), & + J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,EwGIC,E(:),Om(:)) end subroutine GOK_RKS diff --git a/src/eDFT/LIM_RKS.f90 b/src/eDFT/LIM_RKS.f90 index 44da113..7708b1c 100644 --- a/src/eDFT/LIM_RKS.f90 +++ b/src/eDFT/LIM_RKS.f90 @@ -1,5 +1,5 @@ -subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c) +subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight, & + maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c) ! Perform restricted Kohn-Sham calculation for ensembles @@ -10,6 +10,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres integer,intent(in) :: x_rung,c_rung character(len=12),intent(in) :: x_DFA,c_DFA + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -65,7 +66,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres write(*,'(A40)') '*************************************************' write(*,*) - call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,nEns,wLIM,nGrid,weight,maxSCF,thresh, & + call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,nGrid,weight,maxSCF,thresh, & max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwZW,EwGICZW,c) !------------------------------------------------------------------------ @@ -82,7 +83,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres write(*,'(A40)') '*************************************************' write(*,*) - call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh, & + call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,maxSCF,thresh, & max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwEW,EwGICEW,c) !------------------------------------------------------------------------ diff --git a/src/eDFT/RMFL20_lda_correlation_energy.f90 b/src/eDFT/RMFL20_lda_correlation_energy.f90 index 4de6375..8826761 100644 --- a/src/eDFT/RMFL20_lda_correlation_energy.f90 +++ b/src/eDFT/RMFL20_lda_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec) +subroutine RMFL20_lda_correlation_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ec) ! Compute the restricted version of the Marut-Fromager-Loos weight-dependent correlation functional ! The RMFL20 is a two-state, single-weight correlation functional for spin-unpolarized systems @@ -8,6 +8,7 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec) ! Input variables + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -16,7 +17,6 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec) ! Local variables - logical :: LDA_centered = .true. integer :: iEns double precision :: EcLDA double precision,allocatable :: aMFL(:,:) @@ -55,7 +55,7 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec) EceLDA(2) = EcLDA + wEns(2)*(EceLDA(2) - EceLDA(1)) EceLDA(1) = EcLDA - end if + end if ! Weight-denpendent functional for ensembles diff --git a/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 b/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 index 1423465..f405063 100644 --- a/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 +++ b/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ec) +subroutine RMFL20_lda_correlation_individual_energy(LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,Ec) ! Compute eLDA correlation energy @@ -7,6 +7,7 @@ subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow, ! Input variables + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -16,7 +17,6 @@ subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow, ! Local variables - logical :: LDA_centered = .true. integer :: iEns double precision :: EcLDA double precision,allocatable :: aMFL(:,:) @@ -48,7 +48,6 @@ subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow, ! LDA-centered functional - if(LDA_centered) then call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA) diff --git a/src/eDFT/RMFL20_lda_correlation_potential.f90 b/src/eDFT/RMFL20_lda_correlation_potential.f90 index 19f3e49..0007e35 100644 --- a/src/eDFT/RMFL20_lda_correlation_potential.f90 +++ b/src/eDFT/RMFL20_lda_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine RMFL20_lda_correlation_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) +subroutine RMFL20_lda_correlation_potential(LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) ! Compute Marut-Fromager-Loos weight-dependent LDA correlation potential @@ -7,6 +7,7 @@ subroutine RMFL20_lda_correlation_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,F ! Input variables + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -17,7 +18,6 @@ subroutine RMFL20_lda_correlation_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,F ! Local variables - logical :: LDA_centered = .true. integer :: iEns double precision,allocatable :: aMFL(:,:) double precision,allocatable :: FcLDA(:,:) diff --git a/src/eDFT/RMFL20_lda_exchange_energy.f90 b/src/eDFT/RMFL20_lda_exchange_energy.f90 index 555c257..5e9cb1e 100644 --- a/src/eDFT/RMFL20_lda_exchange_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) +subroutine RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex) ! Compute the restricted version of the Marut-Fromager-Loos weight-dependent exchange functional ! The RMFL20 is a two-state, single-weight exchange functional @@ -8,6 +8,7 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) ! Input variables + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -16,7 +17,6 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) ! Local variables - logical :: LDA_centered = .true. integer :: iG double precision :: Cxw double precision :: r diff --git a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 index 4feab64..6551b39 100644 --- a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex) +subroutine RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,Ex) ! Compute the restricted version of the Marut-Fromager-Loos 2020 weight-dependent exchange functional @@ -7,6 +7,7 @@ subroutine RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho ! Input variables + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -16,7 +17,6 @@ subroutine RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho ! Local variables - logical :: LDA_centered = .true. integer :: iG double precision :: Cxw double precision :: r,rI @@ -43,9 +43,11 @@ subroutine RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho rI = max(0d0,rho(iG)) if(r > threshold .and. rI > threshold) then + e = Cxw*r**(1d0/3d0) dedr = 1d0/3d0*Cxw*r**(-2d0/3d0) Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r) + endif enddo diff --git a/src/eDFT/RMFL20_lda_exchange_potential.f90 b/src/eDFT/RMFL20_lda_exchange_potential.f90 index 989573b..ce6cfca 100644 --- a/src/eDFT/RMFL20_lda_exchange_potential.f90 +++ b/src/eDFT/RMFL20_lda_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) +subroutine RMFL20_lda_exchange_potential(LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) ! Compute the restricted version of the weight-dependent MFL20 exchange potential @@ -7,6 +7,7 @@ subroutine RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) ! Input variables + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -17,7 +18,6 @@ subroutine RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) ! Local variables - logical :: LDA_centered = .true. integer :: mu,nu,iG double precision :: Cxw double precision :: r,vAO diff --git a/src/eDFT/RS51_lda_exchange_individual_energy.f90 b/src/eDFT/RS51_lda_exchange_individual_energy.f90 index 7ecf521..f2cdfad 100644 --- a/src/eDFT/RS51_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RS51_lda_exchange_individual_energy.f90 @@ -30,10 +30,11 @@ subroutine RS51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex) r = max(0d0,rhow(iG)) rI = max(0d0,rho(iG)) - if(r > threshold .and. rI > threshold) then + if(r > threshold .or. rI > threshold) then e = CxLDA*r**(1d0/3d0) dedr = 1d0/3d0*CxLDA*r**(-2d0/3d0) + Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r) endif diff --git a/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 index d0549ec..8ff4a95 100644 --- a/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 +++ b/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 @@ -19,7 +19,7 @@ subroutine RVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) double precision :: r,rI,rs,x double precision :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p double precision :: dxdrs,dxdx_p,decdx_p - double precision :: drsdra,decdra_p + double precision :: drsdr,decdr_p double precision :: ec_p ! Output variables @@ -42,8 +42,6 @@ subroutine RVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) r = max(0d0,rhow(iG)) rI = max(0d0,rho(iG)) -! spin-up contribution - if(r > threshold .and. rI > threshold) then rs = (4d0*pi*r/3d0)**(-1d0/3d0) @@ -56,17 +54,17 @@ subroutine RVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) ) - drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) - dxdrs = 0.5d0/sqrt(rs) + drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) + dxdrs = 0.5d0/sqrt(rs) dxdx_p = 2d0*x + b_p decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p & - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) ) - decdra_p = drsdra*dxdrs*decdx_p + decdr_p = drsdr*dxdrs*decdx_p - Ec = Ec + weight(iG)*(ec_p*rI + decdra_p*r*rI - decdra_p*r*r) + Ec = Ec + weight(iG)*(ec_p*rI + decdr_p*r*rI - decdr_p*r*r) end if diff --git a/src/eDFT/RVWN5_lda_correlation_potential.f90 b/src/eDFT/RVWN5_lda_correlation_potential.f90 index 85054bb..918f18a 100644 --- a/src/eDFT/RVWN5_lda_correlation_potential.f90 +++ b/src/eDFT/RVWN5_lda_correlation_potential.f90 @@ -20,7 +20,7 @@ subroutine RVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) double precision :: r,rs,x double precision :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p double precision :: dxdrs,dxdx_p,decdx_p - double precision :: drsdra,decdra_p + double precision :: drsdr,decdr_p double precision :: ec_p @@ -57,17 +57,17 @@ subroutine RVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) ) - drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) - dxdrs = 0.5d0/sqrt(rs) + drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) + dxdrs = 0.5d0/sqrt(rs) dxdx_p = 2d0*x + b_p decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p & - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) ) - decdra_p = drsdra*dxdrs*decdx_p + decdr_p = drsdr*dxdrs*decdx_p - Fc(mu,nu) = Fc(mu,nu) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_p + decdra_p*r) + Fc(mu,nu) = Fc(mu,nu) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_p + decdr_p*r) end if diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index fd9e70f..993e8c8 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -32,6 +32,7 @@ program eDFT character(len=7) :: method integer :: x_rung,c_rung character(len=12) :: x_DFA ,c_DFA + logical :: LDA_centered = .true. integer :: SGn double precision :: radial_precision @@ -152,6 +153,8 @@ program eDFT allocate(AO(nBas,nGrid),dAO(ncart,nBas,nGrid)) call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,nGrid,root,AO,dAO) + LDA_centered = .true. + !------------------------------------------------------------------------ ! Compute GOK-RKS energy !------------------------------------------------------------------------ @@ -159,9 +162,9 @@ program eDFT if(method == 'GOK-RKS') then call cpu_time(start_KS) - call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & - nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc, & - Ew,EwGIC,c(:,:)) + call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), & + maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & + S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,EwGIC,c(:,:)) call cpu_time(end_KS) t_KS = end_KS - start_KS @@ -177,9 +180,9 @@ program eDFT if(method == 'LIM-RKS') then call cpu_time(start_KS) - call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & - nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc, & - c(:,:)) + call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), & + maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & + S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:)) call cpu_time(end_KS) t_KS = end_KS - start_KS diff --git a/src/eDFT/exchange_energy.f90 b/src/eDFT/exchange_energy.f90 index 6c0b39e..88e67a8 100644 --- a/src/eDFT/exchange_energy.f90 +++ b/src/eDFT/exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine exchange_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,FxHF,rho,drho,Ex) +subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,P,FxHF,rho,drho,Ex) ! Compute the exchange energy @@ -9,6 +9,7 @@ subroutine exchange_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,FxHF,rho,drho, integer,intent(in) :: rung character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -40,7 +41,7 @@ subroutine exchange_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,FxHF,rho,drho, case(1) - call lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,ExLDA) + call lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,ExLDA) Ex = ExLDA diff --git a/src/eDFT/exchange_individual_energy.f90 b/src/eDFT/exchange_individual_energy.f90 index 25cf221..5698bef 100644 --- a/src/eDFT/exchange_individual_energy.f90 +++ b/src/eDFT/exchange_individual_energy.f90 @@ -1,5 +1,5 @@ -subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas, & - ERI,P,rhow,drhow,rho,drho,Ex) +subroutine exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas, & + ERI,Pw,P,rhow,drhow,rho,drho,Ex) ! Compute the exchange individual energy @@ -10,12 +10,14 @@ subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas, & integer,intent(in) :: rung character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) integer,intent(in) :: nBas double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Pw(nBas,nBas) double precision,intent(in) :: P(nBas,nBas) double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: drhow(ncart,nGrid) @@ -43,7 +45,7 @@ subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas, & case(1) - call lda_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,ExLDA) + call lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,ExLDA) Ex = ExLDA @@ -65,7 +67,7 @@ subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas, & case(666) - call fock_exchange_individual_energy(nBas,P,ERI,ExHF) + call fock_exchange_individual_energy(nBas,Pw,P,ERI,ExHF) Ex = ExHF diff --git a/src/eDFT/exchange_potential.f90 b/src/eDFT/exchange_potential.f90 index 4048d7a..758be66 100644 --- a/src/eDFT/exchange_potential.f90 +++ b/src/eDFT/exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine exchange_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,ERI,AO,dAO,rho,drho,Fx,FxHF) +subroutine exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,P,ERI,AO,dAO,rho,drho,Fx,FxHF) ! Compute the exchange potential @@ -9,6 +9,7 @@ subroutine exchange_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,ERI,AO,dAO, integer,intent(in) :: rung character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -44,7 +45,7 @@ subroutine exchange_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,ERI,AO,dAO, case(1) - call lda_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + call lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) ! GGA functionals diff --git a/src/eDFT/fock_exchange_individual_energy.f90 b/src/eDFT/fock_exchange_individual_energy.f90 index 7069d71..0dcb816 100644 --- a/src/eDFT/fock_exchange_individual_energy.f90 +++ b/src/eDFT/fock_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine fock_exchange_individual_energy(nBas,P,ERI,Ex) +subroutine fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex) ! Compute the Fock exchange potential @@ -7,6 +7,7 @@ subroutine fock_exchange_individual_energy(nBas,P,ERI,Ex) ! Input variables integer,intent(in) :: nBas + double precision,intent(in) :: Pw(nBas,nBas) double precision,intent(in) :: P(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) diff --git a/src/eDFT/hartree_individual_energy.f90 b/src/eDFT/hartree_individual_energy.f90 new file mode 100644 index 0000000..248c120 --- /dev/null +++ b/src/eDFT/hartree_individual_energy.f90 @@ -0,0 +1,65 @@ +subroutine hartree_individual_energy(rung,nBas,ERI,J,Pw,P,EJ) + +! Compute the exchange individual energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: rung + integer,intent(in) :: nBas + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: J(nBas,nBas) + double precision,intent(in) :: Pw(nBas,nBas) + double precision,intent(in) :: P(nBas,nBas) + +! Local variables + + double precision,external :: trace_matrix + +! Output variables + + double precision,intent(out) :: EJ + + select case (rung) + +! Hartree calculation + + case(0) + + call hartree_coulomb(nBas,P(:,:),ERI(:,:,:,:),J(:,:)) + EJ = 0.5d0*trace_matrix(nBas,matmul(P(:,:),J(:,:))) + +! LDA functionals + + case(1) + + call hartree_coulomb(nBas,Pw(:,:),ERI(:,:,:,:),J(:,:)) + EJ = trace_matrix(nBas,matmul(P(:,:),J(:,:))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),J(:,:))) + +! GGA functionals + + case(2) + + call print_warning('!!! Hartee individual energies NYI for GGAs !!!') + stop + +! Hybrid functionals + + case(4) + + call print_warning('!!! Hartree individual energies NYI for Hybrids !!!') + stop + +! Hartree-Fock calculation + + case(666) + + call hartree_coulomb(nBas,P(:,:),ERI(:,:,:,:),J(:,:)) + EJ = 0.5d0*trace_matrix(nBas,matmul(P(:,:),J(:,:))) + + end select + +end subroutine hartree_individual_energy diff --git a/src/eDFT/lda_exchange_energy.f90 b/src/eDFT/lda_exchange_energy.f90 index f265c91..4a153df 100644 --- a/src/eDFT/lda_exchange_energy.f90 +++ b/src/eDFT/lda_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,Ex) +subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ex) ! Select LDA exchange functional @@ -8,6 +8,7 @@ subroutine lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,Ex) ! Input variables character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -38,7 +39,7 @@ subroutine lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,Ex) case ('RMFL20') - call RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) + call RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex) case default diff --git a/src/eDFT/lda_exchange_individual_energy.f90 b/src/eDFT/lda_exchange_individual_energy.f90 index 8b80244..0d120f1 100644 --- a/src/eDFT/lda_exchange_individual_energy.f90 +++ b/src/eDFT/lda_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine lda_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,Ex) +subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,Ex) ! Compute LDA exchange energy for individual states @@ -7,6 +7,7 @@ subroutine lda_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,Ex ! Input variables + logical,intent(in) :: LDA_centered character(len=12),intent(in) :: DFA integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) @@ -29,7 +30,7 @@ subroutine lda_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,Ex case ('RMFL20') - call RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) + call RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) case default diff --git a/src/eDFT/lda_exchange_potential.f90 b/src/eDFT/lda_exchange_potential.f90 index 7bc8606..d3e5bc9 100644 --- a/src/eDFT/lda_exchange_potential.f90 +++ b/src/eDFT/lda_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine lda_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) +subroutine lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) ! Select LDA correlation potential @@ -8,6 +8,7 @@ subroutine lda_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) ! Input variables + logical,intent(in) :: LDA_centered character(len=12),intent(in) :: DFA integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) @@ -41,7 +42,7 @@ subroutine lda_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) case ('RMFL20') - call RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + call RMFL20_lda_exchange_potential(LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) case default diff --git a/src/eDFT/restricted_correlation_energy.f90 b/src/eDFT/restricted_correlation_energy.f90 index ec8491b..d44a838 100644 --- a/src/eDFT/restricted_correlation_energy.f90 +++ b/src/eDFT/restricted_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine restricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) +subroutine restricted_correlation_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,drho,Ec) ! Compute the correlation energy @@ -9,6 +9,7 @@ subroutine restricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drh integer,intent(in) :: rung character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -38,7 +39,7 @@ subroutine restricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drh case(1) - call restricted_lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:),Ec) + call restricted_lda_correlation_energy(DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),rho(:),Ec) ! GGA functionals diff --git a/src/eDFT/restricted_correlation_individual_energy.f90 b/src/eDFT/restricted_correlation_individual_energy.f90 index e6cfdb3..93538c8 100644 --- a/src/eDFT/restricted_correlation_individual_energy.f90 +++ b/src/eDFT/restricted_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine restricted_correlation_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ec) +subroutine restricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ec) ! Compute the correlation energy of individual states @@ -9,6 +9,7 @@ subroutine restricted_correlation_individual_energy(rung,DFA,nEns,wEns,nGrid,wei integer,intent(in) :: rung character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -40,7 +41,7 @@ subroutine restricted_correlation_individual_energy(rung,DFA,nEns,wEns,nGrid,wei case(1) - call restricted_lda_correlation_individual_energy(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),rho(:),Ec) + call restricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),rhow(:),rho(:),Ec) ! GGA functionals diff --git a/src/eDFT/restricted_correlation_potential.f90 b/src/eDFT/restricted_correlation_potential.f90 index a731203..78cfbc3 100644 --- a/src/eDFT/restricted_correlation_potential.f90 +++ b/src/eDFT/restricted_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine restricted_correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) +subroutine restricted_correlation_potential(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! Compute the correlation potential @@ -9,6 +9,7 @@ subroutine restricted_correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas integer,intent(in) :: rung character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -41,7 +42,7 @@ subroutine restricted_correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas case(1) - call restricted_lda_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:)) + call restricted_lda_correlation_potential(DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:)) ! GGA functionals diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index f0a9467..ee13f40 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas, & +subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas, & nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,EwGIC,E,Om) ! Compute individual energies as well as excitation energies @@ -10,6 +10,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri integer,intent(in) :: x_rung,c_rung character(len=12),intent(in) :: x_DFA,c_DFA + logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid @@ -75,12 +76,11 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri end do !------------------------------------------------------------------------ -! Hartree energy +! Individua Hartree energy !------------------------------------------------------------------------ do iEns=1,nEns - call hartree_coulomb(nBas,P(:,:,iEns),ERI(:,:,:,:),J(:,:)) - EJ(iEns) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,iEns),J(:,:))) + call hartree_individual_energy(x_rung,nBas,ERI,J(:,:),Pw(:,:),P(:,:,iEns),EJ(iEns)) end do !------------------------------------------------------------------------ @@ -88,10 +88,8 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri !------------------------------------------------------------------------ do iEns=1,nEns - - call exchange_individual_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,ERI(:,:,:,:), & - P(:,:,iEns),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns),Ex(iEns)) - + call exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas,ERI(:,:,:,:), & + Pw(:,:),P(:,:,iEns),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns),Ex(iEns)) end do !------------------------------------------------------------------------ @@ -99,15 +97,10 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri !------------------------------------------------------------------------ do iEns=1,nEns - - call restricted_correlation_individual_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:), & + call restricted_correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:), & rho(:,iEns),drho(:,:,iEns),Ec(iEns)) - end do - - - !------------------------------------------------------------------------ ! Compute auxiliary energies !------------------------------------------------------------------------ diff --git a/src/eDFT/restricted_lda_correlation_energy.f90 b/src/eDFT/restricted_lda_correlation_energy.f90 index 0c7aa9e..7adcfa3 100644 --- a/src/eDFT/restricted_lda_correlation_energy.f90 +++ b/src/eDFT/restricted_lda_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine restricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) +subroutine restricted_lda_correlation_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ec) ! Select LDA correlation functional @@ -7,6 +7,7 @@ subroutine restricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) ! Input variables + logical,intent(in) :: LDA_centered character(len=12),intent(in) :: DFA integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) @@ -34,7 +35,7 @@ subroutine restricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) case ('RMFL20') - call RMFL20_lda_correlation_energy(nEns,wEns(:),nGrid,weight(:),rho(:),Ec) + call RMFL20_lda_correlation_energy(LDA_centered,nEns,wEns(:),nGrid,weight(:),rho(:),Ec) case default diff --git a/src/eDFT/restricted_lda_correlation_individual_energy.f90 b/src/eDFT/restricted_lda_correlation_individual_energy.f90 index e6dd87d..a9c2acd 100644 --- a/src/eDFT/restricted_lda_correlation_individual_energy.f90 +++ b/src/eDFT/restricted_lda_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine restricted_lda_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,Ec) +subroutine restricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,Ec) ! Compute LDA correlation energy for individual states @@ -7,6 +7,7 @@ subroutine restricted_lda_correlation_individual_energy(DFA,nEns,wEns,nGrid,weig ! Input variables + logical,intent(in) :: LDA_centered character(len=12),intent(in) :: DFA integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) @@ -33,7 +34,7 @@ subroutine restricted_lda_correlation_individual_energy(DFA,nEns,wEns,nGrid,weig case ('RMFL20') - call RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ec) + call RMFL20_lda_correlation_individual_energy(LDA_centered,nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ec) case default diff --git a/src/eDFT/restricted_lda_correlation_potential.f90 b/src/eDFT/restricted_lda_correlation_potential.f90 index 0872f2b..108346a 100644 --- a/src/eDFT/restricted_lda_correlation_potential.f90 +++ b/src/eDFT/restricted_lda_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine restricted_lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) +subroutine restricted_lda_correlation_potential(DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) ! Select LDA correlation potential @@ -7,6 +7,7 @@ subroutine restricted_lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas, ! Input variables + logical,intent(in) :: LDA_centered character(len=12),intent(in) :: DFA integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) @@ -36,7 +37,7 @@ subroutine restricted_lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas, case ('RMFL20') - call RMFL20_lda_correlation_potential(nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:)) + call RMFL20_lda_correlation_potential(LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:)) case default