DD eLDA

input/dft

@@ -1,14 +1,24 @@
 # Restricted or unrestricted KS calculation
   GOK-RKS
-# exchange rung: Hartree = 0, LDA = 1 (RS51,S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666
+# exchange rung: 
-  1 RS51
+# Hartree      = 0
-# correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666
+# LDA          = 1: RS51,S51,RMFL20
+# GGA          = 2: G96,B88
+# Hybrid       = 4 
+# Hartree-Fock = 666
+  1 RMFL20
+# correlation rung: 
+# Hartree      = 0
+# LDA          = 1: W38,VWN5,C16,RMFL20
+# GGA          = 2: LYP
+# Hybrid       = 4: B3LYP
+# Hartree-Fock = 666
   1 RVWN5
 # quadrature grid SG-n
   1
 # Number of states in ensemble (nEns)
   2
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
-  0.00000 0.00000
+  0.50000 0.00000
-# eKS: maxSCF thresh   DIIS n_diis guess_type ortho_type
+# GOK-DFT: maxSCF thresh   DIIS n_diis guess_type ortho_type
            64    0.0000001  T     5      1          1

src/eDFT/RMFL20_lda_exchange_energy.f90

@@ -33,7 +33,7 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
   Cx1   = - (176d0/105d0)*(1d0/pi)**(1d0/3d0)
   CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0)
 
-  Cxw   = CxLDA + wEns(1)*(Cx1 - Cx0)
+  Cxw   = CxLDA + wEns(2)*(Cx1 - Cx0)
 
 ! Compute LDA exchange energy
 

src/eDFT/RMFL20_lda_exchange_potential.f90

@@ -28,13 +28,13 @@ subroutine RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
 
   double precision,intent(out)  :: Fx(nBas,nBas)
 
-! Cx coefficient for Slater LDA exchange
+! Weight-dependent Cx coefficient for RMFL20 exchange functional
 
   Cx0   = -(4d0/3d0)*(1d0/pi)**(1d0/3d0)
   Cx1   = -(176d0/105d0)*(1d0/pi)**(1d0/3d0)
   CxLDA = -(3d0/4d0)*(3d0/pi)**(1d0/3d0)
 
-  Cxw   = CxLDA + wEns(1)*(Cx1 - Cx0)
+  Cxw   = CxLDA + wEns(2)*(Cx1 - Cx0)
 
 ! Compute LDA exchange matrix in the AO basis
 

src/eDFT/lda_exchange_derivative_discontinuity.f90

@@ -35,8 +35,7 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow
 
     case ('RMFL20')
 
-!     call MFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:))
+      call RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:))
-      ExDD(:) = 0d0
 
     case default
 

src/eDFT/lda_exchange_individual_energy.f90

@@ -29,7 +29,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(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex)
 
     case default
 

src/eDFT/print_restricted_individual_energy.f90

@@ -1,4 +1,5 @@
-subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,ExcDD,E,Om)
+subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,ExcDD,E, & 
+                                              Om,Omx,Omc,Omxc,OmxDD,OmcDD,OmxcDD)
 
 ! Print individual energies for eDFT calculation
 
@@ -13,6 +14,8 @@ subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,
   double precision,intent(in)        :: EJ(nEns)
   double precision,intent(in)        :: Ex(nEns),   Ec(nEns),   Exc(nEns)
   double precision,intent(in)        :: ExDD(nEns), EcDD(nEns), ExcDD(nEns)
+  double precision,intent(in)        :: Omx(nEns),  Omc(nEns),  Omxc(nEns)
+  double precision,intent(in)        :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns)
   double precision,intent(in)        :: E(nEns)
   double precision,intent(in)        :: Om(nEns)
 
@@ -25,7 +28,7 @@ subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,
 !------------------------------------------------------------------------
 
   write(*,'(A60)')           '-------------------------------------------------'
-  write(*,'(A50)')           ' Individual Kinetic     energies'
+  write(*,'(A50)')           ' Individual kinetic     energies'
   write(*,'(A60)')           '-------------------------------------------------'
   do iEns=1,nEns
     write(*,'(A40,I2,A2,F16.10,A3)') ' Kinetic     energy state ',iEns,': ',ET(iEns),' au'
@@ -38,7 +41,7 @@ subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,
 !------------------------------------------------------------------------
 
   write(*,'(A60)')           '-------------------------------------------------'
-  write(*,'(A50)')           ' Individual Potential   energies'
+  write(*,'(A50)')           ' Individual potential   energies'
   write(*,'(A60)')           '-------------------------------------------------'
   do iEns=1,nEns
     write(*,'(A40,I2,A2,F16.10,A3)') ' Potential   energy state ',iEns,': ',EV(iEns),' au'
@@ -93,12 +96,12 @@ subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,
   write(*,'(A50)')           ' Derivative discontinuities (DD)    '
   write(*,'(A60)')           '-------------------------------------------------'
   do iEns=1,nEns
-    write(*,'(A40,I2,A2,F16.10,A3)') '  x ensemble derivative ',iEns,': ',ExDD(iEns), ' au'
+    write(*,'(A40,I2,A2,F16.10,A3)') '  x ensemble derivative state ',iEns,': ',ExDD(iEns), ' au'
-    write(*,'(A40,I2,A2,F16.10,A3)') '  c ensemble derivative ',iEns,': ',EcDD(iEns), ' au'
+    write(*,'(A40,I2,A2,F16.10,A3)') '  c ensemble derivative state ',iEns,': ',EcDD(iEns), ' au'
-    write(*,'(A40,I2,A2,F16.10,A3)') ' xc ensemble derivative ',iEns,': ',ExcDD(iEns),' au'
+    write(*,'(A40,I2,A2,F16.10,A3)') ' xc ensemble derivative state ',iEns,': ',ExcDD(iEns),' au'
+    write(*,*)
   end do
   write(*,'(A60)')           '-------------------------------------------------'
-  write(*,*)
 
 !------------------------------------------------------------------------
 ! Total and Excitation energies
@@ -111,12 +114,30 @@ subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,
     write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns),' au'
   end do
   write(*,'(A60)')           '-------------------------------------------------'
+
   do iEns=2,nEns
     write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy  1 ->',iEns,': ',Om(iEns),    ' au'
+    write(*,*)
+    write(*,'(A44,      F16.10,A3)') '     x  energy contribution     : ',Omx(iEns),   ' au'
+    write(*,'(A44,      F16.10,A3)') '     c  energy contribution     : ',Omc(iEns),   ' au'
+    write(*,'(A44,      F16.10,A3)') '    xc  energy contribution     : ',Omxc(iEns),  ' au'
+    write(*,*)
+    write(*,'(A44,      F16.10,A3)') '     x  ensemble derivative     : ',OmxDD(iEns), ' au'
+    write(*,'(A44,      F16.10,A3)') '     c  ensemble derivative     : ',OmcDD(iEns), ' au'
+    write(*,'(A44,      F16.10,A3)') '    xc  ensemble derivative     : ',OmxcDD(iEns),' au'
   end do
   write(*,'(A60)')           '-------------------------------------------------'
+
   do iEns=2,nEns
     write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy  1 ->',iEns,': ',Om(iEns)*HaToeV,    ' eV'
+    write(*,*)
+    write(*,'(A44,      F16.10,A3)') '     x  energy contribution     : ',Omx(iEns)*HaToeV,   ' eV'
+    write(*,'(A44,      F16.10,A3)') '     c  energy contribution     : ',Omc(iEns)*HaToeV,   ' eV'
+    write(*,'(A44,      F16.10,A3)') '    xc  energy contribution     : ',Omxc(iEns)*HaToeV,  ' eV'
+    write(*,*)
+    write(*,'(A44,      F16.10,A3)') '     x  ensemble derivative     : ',OmxDD(iEns)*HaToeV, ' eV'
+    write(*,'(A44,      F16.10,A3)') '     c  ensemble derivative     : ',OmcDD(iEns)*HaToeV, ' eV'
+    write(*,'(A44,      F16.10,A3)') '    xc  ensemble derivative     : ',OmxcDD(iEns)*HaToeV,' eV'
   end do
   write(*,'(A60)')           '-------------------------------------------------'
   write(*,*)

src/eDFT/read_options.f90

@@ -51,11 +51,21 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,th
 
 ! EXCHANGE: read rung of Jacob's ladder
 
+  read(1,*)
+  read(1,*)
+  read(1,*)
+  read(1,*)
+  read(1,*)
   read(1,*)
   read(1,*) x_rung,x_DFA
 
 ! CORRELATION: read rung of Jacob's ladder
 
+  read(1,*)
+  read(1,*)
+  read(1,*)
+  read(1,*)
+  read(1,*)
   read(1,*)
   read(1,*) c_rung,c_DFA
 

src/eDFT/restricted_individual_energy.f90

@@ -44,6 +44,8 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
   double precision              :: EJ(nEns)
   double precision              :: Ex(nEns),   Ec(nEns),   Exc(nEns)
   double precision              :: ExDD(nEns), EcDD(nEns), ExcDD(nEns)
+  double precision              :: Omx(nEns),  Omc(nEns),  Omxc(nEns)
+  double precision              :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns)
 
   double precision,external     :: trace_matrix
 
@@ -126,7 +128,17 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
 !------------------------------------------------------------------------
 
   do iEns=1,nEns
+
     Om(iEns)     = E(iEns)     - E(1)
+
+    Omx(iEns)    = Ex(iEns)    - Ex(1)
+    Omc(iEns)    = Ec(iEns)    - Ec(1)
+    Omxc(iEns)   = Exc(iEns)   - Exc(1)
+
+    OmxDD(iEns)  = ExDD(iEns)  - ExDD(1)
+    OmcDD(iEns)  = EcDD(iEns)  - EcDD(1)
+    OmxcDD(iEns) = ExcDD(iEns) - ExcDD(1)
+
   end do
 
 !------------------------------------------------------------------------
@@ -134,6 +146,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
 !------------------------------------------------------------------------
 
   call print_restricted_individual_energy(nEns,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), &
-                                          ExDD(:),EcDD(:),ExcDD(:),E(:),Om(:))
+                                          ExDD(:),EcDD(:),ExcDD(:),E(:), & 
+                                          Om(:),Omx(:),Omc(:),Omxc(:),OmxDD(:),OmcDD(:),OmxcDD(:))
 
 end subroutine restricted_individual_energy

src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90

@@ -30,7 +30,7 @@ subroutine restricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGr
 
       Ec(:) = 0d0
 
-    case ('VWN5')
+    case ('RVWN5')
 
       Ec(:) = 0d0