From 614bf19a9c11085db52a3fda8e88a17c5de78bb5 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 6 Apr 2020 22:27:13 +0200 Subject: [PATCH] B88 for ensemble starting implementation --- input/basis | 65 +++++++----------- input/dft | 4 +- input/molecule | 5 +- input/molecule.xyz | 5 +- input/weight | 65 +++++++----------- src/eDFT/RB88_gga_exchange_energy.f90 | 2 +- .../RB88_gga_exchange_individual_energy.f90 | 61 +++++++++++++++++ src/eDFT/RB88_gga_exchange_potential.f90 | 2 +- ..._lda_exchange_derivative_discontinuity.f90 | 5 +- src/eDFT/RGIC_lda_exchange_energy.f90 | 4 +- .../RGIC_lda_exchange_individual_energy.f90 | 7 +- src/eDFT/RGIC_lda_exchange_potential.f90 | 4 +- .../exchange_derivative_discontinuity.f90 | 3 +- src/eDFT/exchange_individual_energy.f90 | 6 +- .../gga_exchange_derivative_discontinuity.f90 | 40 +++++++++++ src/eDFT/gga_exchange_energy.f90 | 8 +-- src/eDFT/gga_exchange_individual_energy.f90 | 39 +++++++++++ src/eDFT/gga_exchange_potential.f90 | 8 +-- src/eDFT/hartree_individual_energy.f90 | 67 ------------------- .../lda_exchange_derivative_discontinuity.f90 | 2 +- src/eDFT/lda_exchange_individual_energy.f90 | 2 +- src/eDFT/restricted_individual_energy.f90 | 4 +- 22 files changed, 220 insertions(+), 188 deletions(-) create mode 100644 src/eDFT/RB88_gga_exchange_individual_energy.f90 create mode 100644 src/eDFT/gga_exchange_derivative_discontinuity.f90 create mode 100644 src/eDFT/gga_exchange_individual_energy.f90 delete mode 100644 src/eDFT/hartree_individual_energy.f90 diff --git a/input/basis b/input/basis index dc1936c..f19a2d0 100644 --- a/input/basis +++ b/input/basis @@ -1,43 +1,30 @@ -1 9 -S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 +1 6 +S 8 + 1 17880.0000000 0.0007380 + 2 2683.0000000 0.0056770 + 3 611.5000000 0.0288830 + 4 173.5000000 0.1085400 + 5 56.6400000 0.2909070 + 6 20.4200000 0.4483240 + 7 7.8100000 0.2580260 + 8 1.6530000 0.0150630 +S 8 + 1 17880.0000000 -0.0001720 + 2 2683.0000000 -0.0013570 + 3 611.5000000 -0.0067370 + 4 173.5000000 -0.0276630 + 5 56.6400000 -0.0762080 + 6 20.4200000 -0.1752270 + 7 7.8100000 -0.1070380 + 8 1.6530000 0.5670500 S 1 - 1 0.3258000 1.0000000 -S 1 - 1 0.1027000 1.0000000 -S 1 - 1 0.0252600 1.0000000 + 1 0.4869000 1.0000000 +P 3 + 1 28.3900000 0.0460870 + 2 6.2700000 0.2401810 + 3 1.6950000 0.5087440 P 1 - 1 1.4070000 1.0000000 -P 1 - 1 0.3880000 1.0000000 -P 1 - 1 0.1020000 1.0000000 + 1 0.4317000 1.0000000 D 1 - 1 1.0570000 1.0000000 -D 1 - 1 0.2470000 1.0000000 -2 9 -S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 -S 1 - 1 0.3258000 1.0000000 -S 1 - 1 0.1027000 1.0000000 -S 1 - 1 0.0252600 1.0000000 -P 1 - 1 1.4070000 1.0000000 -P 1 - 1 0.3880000 1.0000000 -P 1 - 1 0.1020000 1.0000000 -D 1 - 1 1.0570000 1.0000000 -D 1 - 1 0.2470000 1.0000000 + 1 2.2020000 1.0000000 diff --git a/input/dft b/input/dft index 7679b2d..617161f 100644 --- a/input/dft +++ b/input/dft @@ -6,14 +6,14 @@ # GGA = 2: RB88 # Hybrid = 4 # Hartree-Fock = 666 - 1 RGIC + 2 RB88 # correlation rung: # Hartree = 0 # LDA = 1: RVWN5,RMFL20 # GGA = 2: # Hybrid = 4: # Hartree-Fock = 666 - 0 H + 1 RVWN5 # quadrature grid SG-n 1 # Number of states in ensemble (nEns) diff --git a/input/molecule b/input/molecule index 8076140..edeba31 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,4 @@ # nAt nEla nElb nCore nRyd - 2 1 1 0 0 + 1 5 5 0 0 # Znuc x y z - H 0.0 0.0 0.0 - H 0.0 0.0 1.4 + Ne 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index 6edc99d..1c70680 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 + Ne 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/weight b/input/weight index dc1936c..f19a2d0 100644 --- a/input/weight +++ b/input/weight @@ -1,43 +1,30 @@ -1 9 -S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 +1 6 +S 8 + 1 17880.0000000 0.0007380 + 2 2683.0000000 0.0056770 + 3 611.5000000 0.0288830 + 4 173.5000000 0.1085400 + 5 56.6400000 0.2909070 + 6 20.4200000 0.4483240 + 7 7.8100000 0.2580260 + 8 1.6530000 0.0150630 +S 8 + 1 17880.0000000 -0.0001720 + 2 2683.0000000 -0.0013570 + 3 611.5000000 -0.0067370 + 4 173.5000000 -0.0276630 + 5 56.6400000 -0.0762080 + 6 20.4200000 -0.1752270 + 7 7.8100000 -0.1070380 + 8 1.6530000 0.5670500 S 1 - 1 0.3258000 1.0000000 -S 1 - 1 0.1027000 1.0000000 -S 1 - 1 0.0252600 1.0000000 + 1 0.4869000 1.0000000 +P 3 + 1 28.3900000 0.0460870 + 2 6.2700000 0.2401810 + 3 1.6950000 0.5087440 P 1 - 1 1.4070000 1.0000000 -P 1 - 1 0.3880000 1.0000000 -P 1 - 1 0.1020000 1.0000000 + 1 0.4317000 1.0000000 D 1 - 1 1.0570000 1.0000000 -D 1 - 1 0.2470000 1.0000000 -2 9 -S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 -S 1 - 1 0.3258000 1.0000000 -S 1 - 1 0.1027000 1.0000000 -S 1 - 1 0.0252600 1.0000000 -P 1 - 1 1.4070000 1.0000000 -P 1 - 1 0.3880000 1.0000000 -P 1 - 1 0.1020000 1.0000000 -D 1 - 1 1.0570000 1.0000000 -D 1 - 1 0.2470000 1.0000000 + 1 2.2020000 1.0000000 diff --git a/src/eDFT/RB88_gga_exchange_energy.f90 b/src/eDFT/RB88_gga_exchange_energy.f90 index 7f7d48a..04c87f6 100644 --- a/src/eDFT/RB88_gga_exchange_energy.f90 +++ b/src/eDFT/RB88_gga_exchange_energy.f90 @@ -38,7 +38,7 @@ subroutine RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) r = max(0d0,0.5d0*rho(iG)) if(r > threshold) then - g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 + g = 0.25d0*(drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2) x = sqrt(g)/r**(4d0/3d0) Ex = Ex + weight(iG)*alpha*r**(4d0/3d0) & diff --git a/src/eDFT/RB88_gga_exchange_individual_energy.f90 b/src/eDFT/RB88_gga_exchange_individual_energy.f90 new file mode 100644 index 0000000..761c13e --- /dev/null +++ b/src/eDFT/RB88_gga_exchange_individual_energy.f90 @@ -0,0 +1,61 @@ +subroutine RB88_gga_exchange_individual_energy(nGrid,weight,rhow,drhow,rho,drho,Ex) + +! Compute restricted Becke's GGA indivudal energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: drhow(ncart,nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(ncart,nGrid) + +! Local variables + + integer :: iG + double precision :: alpha + double precision :: beta + double precision :: r,rI,g,x + double precision :: ex_p,dexdr_p + +! Output variables + + double precision,intent(out) :: Ex + +! Coefficients for B88 GGA exchange functional + + alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) + beta = 0.0042d0 + +! Compute GGA exchange matrix in the AO basis + + Ex = 0d0 + + do iG=1,nGrid + + r = max(0d0,0.5d0*rhow(iG)) + rI = max(0d0,0.5d0*rho(iG)) + + if(r > threshold .and. rI > threshold) then + + g = 0.25d0*(drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2) + x = sqrt(g)/r**(4d0/3d0) + + dexdr_p = 4d0/3d0*r**(1d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2) & + + 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0) & + - 2d0*3d0/4d0*beta*g**(-1d0/4d0)/r**(2d0/3d0) + + ex_p = alpha*r**(4d0/3d0) & + - weight(iG)*beta*x**2*r**(4d0/3d0)/(1d0 + 6d0*beta*x*asinh(x)) + + Ex = Ex + weight(iG)*(ex_p*rI + dexdr_p*r*rI - dexdr_p*r*r) + + end if + + end do + +end subroutine RB88_gga_exchange_individual_energy diff --git a/src/eDFT/RB88_gga_exchange_potential.f90 b/src/eDFT/RB88_gga_exchange_potential.f90 index 48466f6..1dca23b 100644 --- a/src/eDFT/RB88_gga_exchange_potential.f90 +++ b/src/eDFT/RB88_gga_exchange_potential.f90 @@ -43,7 +43,7 @@ subroutine RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) if(r > threshold) then - g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 + g = 0.25d0*(drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2) vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) Fx(mu,nu) = Fx(mu,nu) & + vAO*(4d0/3d0*r**(1d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2) & diff --git a/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 index cd0fa3f..3f27dd2 100644 --- a/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 @@ -35,7 +35,8 @@ subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rho c = - 0.36718902716347124d0 w = wEns(2) - dCxGICdw = CxLDA*(0.5d0*b + (2d0*a + 0.5d0*c)*(w - 0.5d0) - (1d0 - w)*w*(3d0*b + 4d0*c*(w - 0.5d0))) + dCxGICdw = (0.5d0*b + (2d0*a + 0.5d0*c)*(w - 0.5d0) - (1d0 - w)*w*(3d0*b + 4d0*c*(w - 0.5d0))) + dCxGICdw = CxLDA*dCxGICdw dExdw(:) = 0d0 @@ -52,8 +53,6 @@ subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rho end do - ExDD(:) = 0d0 - do iEns=1,nEns do jEns=2,nEns diff --git a/src/eDFT/RGIC_lda_exchange_energy.f90 b/src/eDFT/RGIC_lda_exchange_energy.f90 index dafc17b..f8f0460 100644 --- a/src/eDFT/RGIC_lda_exchange_energy.f90 +++ b/src/eDFT/RGIC_lda_exchange_energy.f90 @@ -32,7 +32,8 @@ subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) c = - 0.36718902716347124d0 w = wEns(2) - CxGIC = CxLDA*w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) + CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) + CxGIC = CxLDA*CxGIC ! Compute GIC-LDA exchange energy @@ -43,7 +44,6 @@ subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) r = max(0d0,rho(iG)) if(r > threshold) then - Ex = Ex + weight(iG)*CxLDA*r**(4d0/3d0) Ex = Ex + weight(iG)*CxGIC*r**(4d0/3d0) endif diff --git a/src/eDFT/RGIC_lda_exchange_individual_energy.f90 b/src/eDFT/RGIC_lda_exchange_individual_energy.f90 index 68c553a..c73c9bc 100644 --- a/src/eDFT/RGIC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RGIC_lda_exchange_individual_energy.f90 @@ -34,7 +34,8 @@ subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,E c = - 0.36718902716347124d0 w = wEns(2) - CxGIC = CxLDA*w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) + CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) + CxGIC = CxLDA*CxGIC ! Compute LDA exchange matrix in the AO basis @@ -46,10 +47,6 @@ subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,E if(r > threshold .and. rI > threshold) then - e_p = CxLDA*r**(1d0/3d0) - dedr = 1d0/3d0*CxLDA*r**(-2d0/3d0) - Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r) - e_p = CxGIC*r**(1d0/3d0) dedr = 1d0/3d0*CxGIC*r**(-2d0/3d0) Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r) diff --git a/src/eDFT/RGIC_lda_exchange_potential.f90 b/src/eDFT/RGIC_lda_exchange_potential.f90 index 0b99764..5ed81a1 100644 --- a/src/eDFT/RGIC_lda_exchange_potential.f90 +++ b/src/eDFT/RGIC_lda_exchange_potential.f90 @@ -34,7 +34,8 @@ subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) c = - 0.36718902716347124d0 w = wEns(2) - CxGIC = CxLDA*w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) + CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) + CxGIC = CxLDA*CxGIC ! Compute LDA exchange matrix in the AO basis @@ -49,7 +50,6 @@ subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) if(r > threshold) then vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) - Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*CxLDA*r**(1d0/3d0) Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*CxGIC*r**(1d0/3d0) endif diff --git a/src/eDFT/exchange_derivative_discontinuity.f90 b/src/eDFT/exchange_derivative_discontinuity.f90 index 223d9e1..088e616 100644 --- a/src/eDFT/exchange_derivative_discontinuity.f90 +++ b/src/eDFT/exchange_derivative_discontinuity.f90 @@ -41,8 +41,7 @@ subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rho case(2) - call print_warning('!!! exchange part of the derivative discontinuity NYI for GGAs !!!') - stop + call gga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:)) ! Hybrid functionals diff --git a/src/eDFT/exchange_individual_energy.f90 b/src/eDFT/exchange_individual_energy.f90 index 5698bef..ee81b97 100644 --- a/src/eDFT/exchange_individual_energy.f90 +++ b/src/eDFT/exchange_individual_energy.f90 @@ -27,6 +27,7 @@ subroutine exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weig ! Local variables double precision :: ExLDA + double precision :: ExGGA double precision :: ExHF ! Output variables @@ -53,8 +54,9 @@ subroutine exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weig case(2) - call print_warning('!!! Individual energies NYI for GGAs !!!') - stop + call gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,ExGGA) + + Ex = ExGGA ! Hybrid functionals diff --git a/src/eDFT/gga_exchange_derivative_discontinuity.f90 b/src/eDFT/gga_exchange_derivative_discontinuity.f90 new file mode 100644 index 0000000..eeedc32 --- /dev/null +++ b/src/eDFT/gga_exchange_derivative_discontinuity.f90 @@ -0,0 +1,40 @@ +subroutine gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD) + +! Compute the exchange GGA part of the derivative discontinuity + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: drhow(ncart,nGrid) + +! Local variables + + +! Output variables + + double precision,intent(out) :: ExDD(nEns) + +! Select correlation functional + + select case (DFA) + + case ('RB88') + + ExDD(:) = 0d0 + + case default + + call print_warning('!!! GGA exchange derivative discontinuity not available !!!') + stop + + end select + +end subroutine gga_exchange_derivative_discontinuity diff --git a/src/eDFT/gga_exchange_energy.f90 b/src/eDFT/gga_exchange_energy.f90 index b110a5b..e81a1af 100644 --- a/src/eDFT/gga_exchange_energy.f90 +++ b/src/eDFT/gga_exchange_energy.f90 @@ -22,27 +22,21 @@ subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) select case (DFA) -! Gill's 96 exchange functional - case ('G96') call G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) -! Becke's 88 exchange functional - case ('RB88') call RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) -! Becke's 88 exchange functional - case ('B88') call B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) case default - call print_warning('!!! GGA exchange functional not available !!!') + call print_warning('!!! GGA exchange energy not available !!!') stop end select diff --git a/src/eDFT/gga_exchange_individual_energy.f90 b/src/eDFT/gga_exchange_individual_energy.f90 new file mode 100644 index 0000000..d233066 --- /dev/null +++ b/src/eDFT/gga_exchange_individual_energy.f90 @@ -0,0 +1,39 @@ +subroutine gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ex) + +! Compute GGA exchange energy for individual states + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: drhow(ncart,nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(ncart,nGrid) + +! Output variables + + double precision :: Ex + +! Select correlation functional + + select case (DFA) + + case ('RB88') + + call RB88_gga_exchange_individual_energy(nGrid,weight(:),rhow(:),drhow(:,:),rho(:),drho(:,:),Ex) + + case default + + call print_warning('!!! GGA exchange individual energy not available !!!') + stop + + end select + +end subroutine gga_exchange_individual_energy diff --git a/src/eDFT/gga_exchange_potential.f90 b/src/eDFT/gga_exchange_potential.f90 index 78f94a2..2e9867d 100644 --- a/src/eDFT/gga_exchange_potential.f90 +++ b/src/eDFT/gga_exchange_potential.f90 @@ -26,27 +26,21 @@ subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drh select case (DFA) -! Gill's 96 exchange functional - case ('G96') call G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) -! Becke's 88 exchange functional - case ('RB88') call RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) -! Becke's 88 exchange functional - case ('B88') call B88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) case default - call print_warning('!!! GGA exchange functional not available !!!') + call print_warning('!!! GGA exchange potential not available !!!') stop end select diff --git a/src/eDFT/hartree_individual_energy.f90 b/src/eDFT/hartree_individual_energy.f90 deleted file mode 100644 index cf1dd3f..0000000 --- a/src/eDFT/hartree_individual_energy.f90 +++ /dev/null @@ -1,67 +0,0 @@ -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,Pw(:,:),ERI(:,:,:,:),J(:,:)) - EJ = trace_matrix(nBas,matmul(P(:,:),J(:,:))) & - - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),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,Pw(:,:),ERI(:,:,:,:),J(:,:)) - EJ = trace_matrix(nBas,matmul(P(:,:),J(:,:))) & - - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),J(:,:))) - - end select - -end subroutine hartree_individual_energy diff --git a/src/eDFT/lda_exchange_derivative_discontinuity.f90 b/src/eDFT/lda_exchange_derivative_discontinuity.f90 index 618b6ba..0f0dac6 100644 --- a/src/eDFT/lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/lda_exchange_derivative_discontinuity.f90 @@ -43,7 +43,7 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow case default - call print_warning('!!! LDA exchange functional not available !!!') + call print_warning('!!! LDA exchange derivative discontinuity not available !!!') stop end select diff --git a/src/eDFT/lda_exchange_individual_energy.f90 b/src/eDFT/lda_exchange_individual_energy.f90 index 1a90522..624c6eb 100644 --- a/src/eDFT/lda_exchange_individual_energy.f90 +++ b/src/eDFT/lda_exchange_individual_energy.f90 @@ -38,7 +38,7 @@ subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weigh case default - call print_warning('!!! LDA correlation functional not available !!!') + call print_warning('!!! LDA exchange individual energy not available !!!') stop end select diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index ee13f40..e501ce3 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -80,7 +80,9 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n !------------------------------------------------------------------------ do iEns=1,nEns - call hartree_individual_energy(x_rung,nBas,ERI,J(:,:),Pw(:,:),P(:,:,iEns),EJ(iEns)) + call hartree_coulomb(nBas,Pw(:,:),ERI(:,:,:,:),J(:,:)) + EJ(iEns) = trace_matrix(nBas,matmul(P(:,:,iEns),J(:,:))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),J(:,:))) end do !------------------------------------------------------------------------