Three-state extension and MOM

examples/molecule.H2

@@ -2,4 +2,4 @@
  2   1     1   0     0
 # Znuc   x            y           z
   H     0.           0.         0.
-  H     0.           0.         2.5
+  H     0.           0.         1.4

input/basis

@@ -1,66 +1,18 @@
-1   21 
+1   3
-S   10
+S   3
-  1  54620.0000000              0.0000180        
+  1     13.0100000              0.0196850        
-  2   8180.0000000              0.0001380        
+  2      1.9620000              0.1379770        
-  3   1862.0000000              0.0007230        
+  3      0.4446000              0.4781480        
-  4    527.3000000              0.0030390        
-  5    172.0000000              0.0109080        
-  6     62.1000000              0.0340350        
-  7     24.2100000              0.0911930        
-  8      9.9930000              0.1992680        
-  9      4.3050000              0.3293550        
- 10      1.9210000              0.3404890        
-S   10
-  1  54620.0000000             -0.0000030        
-  2   8180.0000000             -0.0000250        
-  3   1862.0000000             -0.0001310        
-  4    527.3000000             -0.0005580        
-  5    172.0000000             -0.0019880        
-  6     62.1000000             -0.0063700        
-  7     24.2100000             -0.0172170        
-  8      9.9930000             -0.0408580        
-  9      4.3050000             -0.0742370        
- 10      1.9210000             -0.1192340        
 S   1
-  1      0.8663000              1.0000000        
+  1      0.1220000              1.0000000        
+P   1
+  1      0.7270000              1.0000000
+2   3
+S   3
+  1     13.0100000              0.0196850        
+  2      1.9620000              0.1379770        
+  3      0.4446000              0.4781480        
 S   1
-  1      0.2475000              1.0000000        
+  1      0.1220000              1.0000000        
-S   1
-  1      0.1009000              1.0000000        
-S   1
-  1      0.0412900              1.0000000        
-P   4
-  1     43.7500000              0.0006330        
-  2     10.3300000              0.0048080        
-  3      3.2260000              0.0205270        
-  4      1.1270000              0.0678160        
 P   1
-  1      0.4334000              1.0000000        
+  1      0.7270000              1.0000000
-P   1
-  1      0.1808000              1.0000000        
-P   1
-  1      0.0782700              1.0000000        
-P   1
-  1      0.0337200              1.0000000        
-D   1
-  1      1.6350000              1.0000000        
-D   1
-  1      0.7410000              1.0000000        
-D   1
-  1      0.3350000              1.0000000        
-D   1
-  1      0.1519000              1.0000000        
-F   1
-  1      0.6860000              1.0000000        
-F   1
-  1      0.4010000              1.0000000        
-F   1
-  1      0.2350000              1.0000000        
-G   1
-  1      0.6030000              1.0000000        
-G   1
-  1      0.3240000              1.0000000        
-H   1
-  1      0.5100000              1.0000000
-
-

input/dft

@@ -1,5 +1,5 @@
 # Restricted or unrestricted KS calculation
-  GOK-RKS
+  MOM-RKS
 # exchange rung:
 # Hartree      = 0
 # LDA          = 1: RS51,RMFL20
@@ -13,12 +13,12 @@
 # GGA          = 2: 
 # Hybrid       = 4: 
 # Hartree-Fock = 666
-  1 RMFL20
+  1 RVWN5
 # quadrature grid SG-n
   1
 # Number of states in ensemble (nEns)
-  2
+  3
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
-  0.0
+  0.0000000 0.0000000
 # GOK-DFT: maxSCF thresh   DIIS n_diis guess_type ortho_type
            32    0.00001   T     5      1          1

input/molecule

@@ -1,4 +1,5 @@
 # nAt nEla nElb nCore nRyd
-  1   2     2   0     0
+ 2   1     1   0     0
 # Znuc   x            y           z
-  Be      0.0          0.0         0.0
+  H     0.           0.         0.
+  H     0.           0.         1.4

input/molecule.xyz

@@ -1,3 +1,4 @@
-  1
+  2
 
- Be     0.0000000000    0.0000000000    0.0000000000
+ H      0.0000000000    0.0000000000    0.0000000000
+ H      0.0000000000    0.0000000000    0.7408481486

input/weight

@@ -1,66 +1,18 @@
-1   21 
+1   3
-S   10
+S   3
-  1  54620.0000000              0.0000180        
+  1     13.0100000              0.0196850        
-  2   8180.0000000              0.0001380        
+  2      1.9620000              0.1379770        
-  3   1862.0000000              0.0007230        
+  3      0.4446000              0.4781480        
-  4    527.3000000              0.0030390        
-  5    172.0000000              0.0109080        
-  6     62.1000000              0.0340350        
-  7     24.2100000              0.0911930        
-  8      9.9930000              0.1992680        
-  9      4.3050000              0.3293550        
- 10      1.9210000              0.3404890        
-S   10
-  1  54620.0000000             -0.0000030        
-  2   8180.0000000             -0.0000250        
-  3   1862.0000000             -0.0001310        
-  4    527.3000000             -0.0005580        
-  5    172.0000000             -0.0019880        
-  6     62.1000000             -0.0063700        
-  7     24.2100000             -0.0172170        
-  8      9.9930000             -0.0408580        
-  9      4.3050000             -0.0742370        
- 10      1.9210000             -0.1192340        
 S   1
-  1      0.8663000              1.0000000        
+  1      0.1220000              1.0000000        
+P   1
+  1      0.7270000              1.0000000
+2   3
+S   3
+  1     13.0100000              0.0196850        
+  2      1.9620000              0.1379770        
+  3      0.4446000              0.4781480        
 S   1
-  1      0.2475000              1.0000000        
+  1      0.1220000              1.0000000        
-S   1
-  1      0.1009000              1.0000000        
-S   1
-  1      0.0412900              1.0000000        
-P   4
-  1     43.7500000              0.0006330        
-  2     10.3300000              0.0048080        
-  3      3.2260000              0.0205270        
-  4      1.1270000              0.0678160        
 P   1
-  1      0.4334000              1.0000000        
+  1      0.7270000              1.0000000
-P   1
-  1      0.1808000              1.0000000        
-P   1
-  1      0.0782700              1.0000000        
-P   1
-  1      0.0337200              1.0000000        
-D   1
-  1      1.6350000              1.0000000        
-D   1
-  1      0.7410000              1.0000000        
-D   1
-  1      0.3350000              1.0000000        
-D   1
-  1      0.1519000              1.0000000        
-F   1
-  1      0.6860000              1.0000000        
-F   1
-  1      0.4010000              1.0000000        
-F   1
-  1      0.2350000              1.0000000        
-G   1
-  1      0.6030000              1.0000000        
-G   1
-  1      0.3240000              1.0000000        
-H   1
-  1      0.5100000              1.0000000
-
-

src/eDFT/LIM_RKS.f90

@@ -1,4 +1,4 @@
-subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight, &
+subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,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
@@ -12,7 +12,6 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight
   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
@@ -35,8 +34,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight
 ! Local variables
 
   integer                       :: iEns
-  double precision              :: EwZW
+  double precision              :: Ew(nEnS)
-  double precision              :: EwEW
   double precision              :: wLIM(nEns)
   double precision              :: Om(nEns)
 
@@ -49,6 +47,11 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight
   write(*,*)'************************************************'
   write(*,*)
 
+! Initializatio
+
+  Ew(:) = 0d0
+  Om(:) = 0d0
+
 !------------------------------------------------------------------------
 ! Zero-weight calculation
 !------------------------------------------------------------------------
@@ -57,8 +60,9 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight
   write(*,'(A40)')           '            ZERO-WEIGHT CALCULATION              '
   write(*,'(A40)')           '*************************************************'
 
-  wLIM(1)      = 1d0
+  wLIM(1) = 1d0
-  wLIM(2:nEns) = 0d0
+  wLIM(2) = 0d0
+  wLIM(3) = 0d0
 
   do iEns=1,nEns
     write(*,'(A20,I2,A2,F16.10)') ' Weight of state  ',iEns,': ',wLIM(iEns)
@@ -67,43 +71,71 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight
   write(*,*)
 
   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,c)
+               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c)
 
 !------------------------------------------------------------------------
 ! Equiensemble calculation
 !------------------------------------------------------------------------
 
   write(*,'(A40)')           '*************************************************'
-  write(*,'(A40)')           '        NON-ZERO-WEIGHT CALCULATION              '
+  write(*,'(A40)')           '  TWO-STATE   EQUI-WEIGHT CALCULATION            '
   write(*,'(A40)')           '*************************************************'
 
+  wLIM(1) = 0.5d0
+  wLIM(2) = 0.5d0
+  wLIM(3) = 0.0d0
+
   do iEns=1,nEns
-    write(*,'(A20,I2,A2,F16.10)') ' Weight of state  ',iEns,': ',wEns(iEns)
+    write(*,'(A20,I2,A2,F16.10)') ' Weight of state  ',iEns,': ',wLIM(iEns)
   end do
   write(*,'(A40)')           '*************************************************'
   write(*,*)
 
-  call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,maxSCF,thresh, &
+  call GOK_RKS(.true.,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,EwEW,c)
+               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c)
+
+!------------------------------------------------------------------------
+! Equiensemble calculation
+!------------------------------------------------------------------------
+
+  write(*,'(A40)')           '*************************************************'
+  write(*,'(A40)')           '  THREE-STATE EQUI-WEIGHT CALCULATION            '
+  write(*,'(A40)')           '*************************************************'
+
+  wLIM(1) = 1d0/3d0
+  wLIM(2) = 1d0/3d0
+  wLIM(3) = 1d0/3d0
+
+  do iEns=1,nEns
+    write(*,'(A20,I2,A2,F16.10)') ' Weight of state  ',iEns,': ',wLIM(iEns)
+  end do
+  write(*,'(A40)')           '*************************************************'
+  write(*,*)
+
+  call GOK_RKS(.true.,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,Ew(3),c)
+
 
 !------------------------------------------------------------------------
 ! LIM excitation energies
 !------------------------------------------------------------------------
 
-  Om(:)    = 0d0
+  Om(2) = (Ew(2) - Ew(1))/2d0
-  if(wEns(2) > 10d-3) then
+  Om(3) = (Ew(3) - Ew(1))/3d0 + 0.5d0*Om(2)
-    Om(2)    = (EwEW - EwZW)/wEns(2)
-  end if
 
   write(*,'(A60)')           '-------------------------------------------------'
   write(*,'(A60)')           ' LINEAR INTERPOLATION METHOD EXCITATION ENERGIES '
   write(*,'(A60)')           '-------------------------------------------------'
 
-  write(*,'(A44,F16.10,A3)') ' Zero-weight     ensemble energy',EwZW,   ' au'
+  write(*,'(A44,F16.10,A3)') '     Ensemble   energy  #1      ',Ew(1),' au'
-  write(*,'(A44,F16.10,A3)') ' Equi-weight     ensemble energy',EwEW,    ' au'
+  write(*,'(A44,F16.10,A3)') '     Ensemble   energy  #2      ',Ew(2),' au'
+  write(*,'(A44,F16.10,A3)') '     Ensemble   energy  #3      ',Ew(3),' au'
   write(*,'(A60)')           '-------------------------------------------------'
   do iEns=2,nEns
     write(*,'(A40,I2,A2,F16.10,A3)') '    Excitation energy  1 ->',iEns,': ',Om(iEns),       ' au'
+  end do
+  write(*,*)
+  do iEns=2,nEns
     write(*,'(A40,I2,A2,F16.10,A3)') '    Excitation energy  1 ->',iEns,': ',Om(iEns)*HaToeV,' eV'
   end do
   write(*,'(A60)')           '-------------------------------------------------'

src/eDFT/MOM_RKS.f90

@@ -0,0 +1,142 @@
+subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,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
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  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
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  integer,intent(in)            :: maxSCF,max_diis,guess_type
+  double precision,intent(in)   :: thresh
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: AO(nBas,nGrid)
+  double precision,intent(in)   :: dAO(ncart,nBas,nGrid)
+
+  integer,intent(in)            :: nO,nV
+  double precision,intent(in)   :: S(nBas,nBas)
+  double precision,intent(in)   :: T(nBas,nBas)
+  double precision,intent(in)   :: V(nBas,nBas)
+  double precision,intent(in)   :: Hc(nBas,nBas) 
+  double precision,intent(in)   :: X(nBas,nBas) 
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ENuc
+
+  double precision,intent(out)  :: c(nBas,nBas) 
+
+! Local variables
+
+  integer                       :: iEns
+  double precision              :: Ew(nEns)
+  double precision              :: wMOM(nEns)
+  double precision              :: Om(nEns)
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'************************************************'
+  write(*,*)'*          Maximum Overlao method              *'
+  write(*,*)'*          for excitation energies             *'
+  write(*,*)'************************************************'
+  write(*,*)
+
+! Initialization
+
+  Ew(:) = 0d0
+  Om(:) = 0d0
+
+!------------------------------------------------------------------------
+! Zero-weight calculation
+!------------------------------------------------------------------------
+
+  write(*,'(A40)')           '*************************************************'
+  write(*,'(A40)')           '           STATE-SPECIFIC CALCULATION #1         '
+  write(*,'(A40)')           '*************************************************'
+
+  wMOM(1) = 1d0
+  wMOM(2) = 0d0
+  wMOM(3) = 0d0
+
+  do iEns=1,nEns
+    write(*,'(A20,I2,A2,F16.10)') ' Weight of state  ',iEns,': ',wMOM(iEns)
+  end do
+  write(*,'(A40)')           '*************************************************'
+  write(*,*)
+
+  call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, &
+               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c)
+
+!------------------------------------------------------------------------
+! Equiensemble calculation
+!------------------------------------------------------------------------
+
+  write(*,'(A40)')           '*************************************************'
+  write(*,'(A40)')           '           STATE-SPECIFIC CALCULATION #2         '
+  write(*,'(A40)')           '*************************************************'
+
+  wMOM(1) = 0d0
+  wMOM(2) = 1d0
+  wMOM(3) = 0d0
+
+  do iEns=1,nEns
+    write(*,'(A20,I2,A2,F16.10)') ' Weight of state  ',iEns,': ',wMOM(iEns)
+  end do
+  write(*,'(A40)')           '*************************************************'
+  write(*,*)
+
+  call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, &
+               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c)
+
+!------------------------------------------------------------------------
+! Equiensemble calculation
+!------------------------------------------------------------------------
+
+  write(*,'(A40)')           '*************************************************'
+  write(*,'(A40)')           '           STATE-SPECIFIC CALCULATION #3         '
+  write(*,'(A40)')           '*************************************************'
+
+  wMOM(1) = 0d0
+  wMOM(2) = 0d0
+  wMOM(3) = 1.0d0
+
+  do iEns=1,nEns
+    write(*,'(A20,I2,A2,F16.10)') ' Weight of state  ',iEns,': ',wMOM(iEns)
+  end do
+  write(*,'(A40)')           '*************************************************'
+  write(*,*)
+
+  call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, &
+               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c)
+
+!------------------------------------------------------------------------
+! MOM excitation energies
+!------------------------------------------------------------------------
+
+  Om(:) = Ew(:) - Ew(1)
+
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,'(A60)')           ' LINEAR INTERPOLATION METHOD EXCITATION ENERGIES '
+  write(*,'(A60)')           '-------------------------------------------------'
+
+  write(*,'(A44,F16.10,A3)') '     Ensemble   energy  #1      ',Ew(1),' au'
+  write(*,'(A44,F16.10,A3)') '     Ensemble   energy  #2      ',Ew(2),' au'
+  write(*,'(A44,F16.10,A3)') '     Ensemble   energy  #3      ',Ew(3),' au'
+  write(*,'(A60)')           '-------------------------------------------------'
+  do iEns=2,nEns
+    write(*,'(A40,I2,A2,F16.10,A3)') '    Excitation energy  1 ->',iEns,': ',Om(iEns),       ' au'
+  end do
+  write(*,*)
+  do iEns=2,nEns
+    write(*,'(A40,I2,A2,F16.10,A3)') '    Excitation energy  1 ->',iEns,': ',Om(iEns)*HaToeV,' eV'
+  end do
+  write(*,'(A60)')           '-------------------------------------------------'
+
+
+end subroutine MOM_RKS

src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90

@@ -37,9 +37,9 @@ subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rho
 
 ! Parameters for H2 at equilibrium
 
-! a = + 0.5751782560799208d0
+  a = + 0.5751782560799208d0
-! b = - 0.021108186591137282d0
+  b = - 0.021108186591137282d0
-! c = - 0.36718902716347124d0
+  c = - 0.36718902716347124d0
 
 ! Parameters for stretch H2
 
@@ -55,11 +55,11 @@ subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rho
 
 ! Parameters for HNO
 
-  a = 0.0061158387543040335d0
+! a = 0.0061158387543040335d0
-  b = -0.00005968703047293955d0
+! b = -0.00005968703047293955d0
-  c = -0.00001692245714408755d0
+! c = -0.00001692245714408755d0
 
-  w = wEns(2)
+  w = 0.5d0*wEns(2) + wEns(3)
   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
 
@@ -72,7 +72,8 @@ subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rho
     if(r > threshold) then
  
       dExdw(1) = 0d0
-      dExdw(2) = dExdw(2) + weight(iG)*dCxGICdw*r**(4d0/3d0)
+      dExdw(2) = dExdw(2) + 0.5d0*weight(iG)*dCxGICdw*r**(4d0/3d0)
+      dExdw(3) = dExdw(3) + 1.0d0*weight(iG)*dCxGICdw*r**(4d0/3d0)
 
     end if
      

src/eDFT/RGIC_lda_exchange_energy.f90

@@ -29,9 +29,9 @@ subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
 
 ! Parameters for H2 at equilibrium
 
-! a = + 0.5751782560799208d0
+  a = + 0.5751782560799208d0
-! b = - 0.021108186591137282d0
+  b = - 0.021108186591137282d0
-! c = - 0.36718902716347124d0
+  c = - 0.36718902716347124d0
 
 ! Parameters for stretch H2
 
@@ -47,11 +47,11 @@ subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
 
 ! Parameters for HNO
 
-  a = 0.0061158387543040335d0
+! a = 0.0061158387543040335d0
-  b = -0.00005968703047293955d0
+! b = -0.00005968703047293955d0
-  c = -0.00001692245714408755d0
+! c = -0.00001692245714408755d0
 
-  w = wEns(2)
+  w = 0.5d0*wEns(2) + wEns(3)
   CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)
   CxGIC = CxLDA*CxGIC
 

src/eDFT/RGIC_lda_exchange_individual_energy.f90

@@ -32,9 +32,9 @@ subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,E
 
 ! Parameters for H2 at equilibrium
 
-! a = + 0.5751782560799208d0
+  a = + 0.5751782560799208d0
-! b = - 0.021108186591137282d0
+  b = - 0.021108186591137282d0
-! c = - 0.36718902716347124d0
+  c = - 0.36718902716347124d0
 
 ! Parameters for stretch H2
 
@@ -51,11 +51,11 @@ subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,E
 
 ! Parameters for HNO
 
-  a = 0.0061158387543040335d0
+! a = 0.0061158387543040335d0
-  b = -0.00005968703047293955d0
+! b = -0.00005968703047293955d0
-  c = -0.00001692245714408755d0
+! c = -0.00001692245714408755d0
 
-  w = wEns(2)
+  w = 0.5d0*wEns(2) + wEns(3)
   CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)
   CxGIC = CxLDA*CxGIC
  

src/eDFT/RGIC_lda_exchange_potential.f90

@@ -32,9 +32,9 @@ subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
 
 ! Parameters for H2 at equilibrium
 
-! a = + 0.5751782560799208d0
+  a = + 0.5751782560799208d0
-! b = - 0.021108186591137282d0
+  b = - 0.021108186591137282d0
-! c = - 0.36718902716347124d0
+  c = - 0.36718902716347124d0
 
 ! Parameters for stretch H2
 
@@ -51,11 +51,11 @@ subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
 
 ! Parameters for HNO
 
-  a = 0.0061158387543040335d0
+! a = 0.0061158387543040335d0
-  b = -0.00005968703047293955d0
+! b = -0.00005968703047293955d0
-  c = -0.00001692245714408755d0
+! c = -0.00001692245714408755d0
 
-  w = wEns(2)
+  w = 0.5d0*wEns(2) + wEns(3)
   CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)
   CxGIC = CxLDA*CxGIC
   

src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90

@@ -34,9 +34,13 @@ subroutine RMFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weigh
   aMFL(2,1) = +0.00540994d0
   aMFL(3,1) = +0.0830766d0
 
-  aMFL(1,2) = -0.0144633d0
+  aMFL(1,2) = -0.0282814d0
-  aMFL(2,2) = -0.0506019d0
+  aMFL(2,2) = +0.00273925d0
-  aMFL(3,2) = +0.0331417d0
+  aMFL(3,2) = +0.0664914d0
+
+  aMFL(1,3) = -0.0144633d0
+  aMFL(2,3) = -0.0506019d0
+  aMFL(3,3) = +0.0331417d0
 
 ! Compute correlation energy for ground, singly-excited and doubly-excited states
 

src/eDFT/RMFL20_lda_correlation_energy.f90

@@ -36,9 +36,13 @@ subroutine RMFL20_lda_correlation_energy(LDA_centered,nEns,wEns,nGrid,weight,rho
   aMFL(2,1) = +0.00540994d0
   aMFL(3,1) = +0.0830766d0
 
-  aMFL(1,2) = -0.0144633d0
+  aMFL(1,2) = -0.0282814d0
-  aMFL(2,2) = -0.0506019d0
+  aMFL(2,2) = +0.00273925d0
-  aMFL(3,2) = +0.0331417d0
+  aMFL(3,2) = +0.0664914d0
+
+  aMFL(1,3) = -0.0144633d0
+  aMFL(2,3) = -0.0506019d0
+  aMFL(3,3) = +0.0331417d0
 
 ! Compute correlation energy for ground and doubly-excited states
 
@@ -52,6 +56,7 @@ subroutine RMFL20_lda_correlation_energy(LDA_centered,nEns,wEns,nGrid,weight,rho
 
     call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),EcLDA)
 
+    EceLDA(3) = EcLDA + wEns(3)*(EceLDA(3) - EceLDA(1))
     EceLDA(2) = EcLDA + wEns(2)*(EceLDA(2) - EceLDA(1))
     EceLDA(1) = EcLDA
 

src/eDFT/RMFL20_lda_correlation_individual_energy.f90

@@ -36,9 +36,13 @@ subroutine RMFL20_lda_correlation_individual_energy(LDA_centered,nEns,wEns,nGrid
   aMFL(2,1) = +0.00540994d0
   aMFL(3,1) = +0.0830766d0
 
-  aMFL(1,2) = -0.0144633d0
+  aMFL(1,2) = -0.0282814d0
-  aMFL(2,2) = -0.0506019d0
+  aMFL(2,2) = +0.00273925d0
-  aMFL(3,2) = +0.0331417d0
+  aMFL(3,2) = +0.0664914d0
+
+  aMFL(1,3) = -0.0144633d0
+  aMFL(2,3) = -0.0506019d0
+  aMFL(3,3) = +0.0331417d0
 
 ! Compute correlation energy for ground- and doubly-excited states
 
@@ -52,6 +56,7 @@ subroutine RMFL20_lda_correlation_individual_energy(LDA_centered,nEns,wEns,nGrid
   
     call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA)
 
+    EceLDA(3) = EcLDA + wEns(3)*(EceLDA(3) - EceLDA(1))
     EceLDA(2) = EcLDA + wEns(2)*(EceLDA(2) - EceLDA(1))
     EceLDA(1) = EcLDA
 

src/eDFT/RMFL20_lda_correlation_potential.f90

@@ -37,9 +37,13 @@ subroutine RMFL20_lda_correlation_potential(LDA_centered,nEns,wEns,nGrid,weight,
   aMFL(2,1) = +0.00540994d0
   aMFL(3,1) = +0.0830766d0
 
-  aMFL(1,2) = -0.0144633d0
+  aMFL(1,2) = -0.0282814d0
-  aMFL(2,2) = -0.0506019d0
+  aMFL(2,2) = +0.00273925d0
-  aMFL(3,2) = +0.0331417d0
+  aMFL(3,2) = +0.0664914d0
+
+  aMFL(1,3) = -0.0144633d0
+  aMFL(2,3) = -0.0506019d0
+  aMFL(3,3) = +0.0331417d0
 
 ! Compute correlation energy for ground- and doubly-excited states
 
@@ -53,6 +57,7 @@ subroutine RMFL20_lda_correlation_potential(LDA_centered,nEns,wEns,nGrid,weight,
 
     call RVWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:),FcLDA(:,:))
 
+    FceLDA(:,:,3) = FcLDA(:,:) + wEns(3)*(FceLDA(:,:,3) - FceLDA(:,:,1))
     FceLDA(:,:,2) = FcLDA(:,:) + wEns(2)*(FceLDA(:,:,2) - FceLDA(:,:,1))
     FceLDA(:,:,1) = FcLDA(:,:)
 

src/eDFT/eDFT.f90

@@ -174,13 +174,13 @@ program eDFT
  end if
 
 !------------------------------------------------------------------------
-! Compute RKS energy
+! Compute LIM excitation energies
 !------------------------------------------------------------------------
 
   if(method == 'LIM-RKS') then
 
     call cpu_time(start_KS)
-    call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),   &
+    call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,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)
@@ -191,6 +191,24 @@ program eDFT
 
   end if
 
+!------------------------------------------------------------------------
+! Compute MOM excitation energies
+!------------------------------------------------------------------------
+
+  if(method == 'MOM-RKS') then
+
+    call cpu_time(start_KS)
+    call MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,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
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM-RKS = ',t_KS,' seconds'
+    write(*,*)
+
+  end if
+
 !------------------------------------------------------------------------
 ! Compute GOK-UKS energy (BROKEN)
 !------------------------------------------------------------------------

src/eDFT/print_restricted_individual_energy.f90

@@ -147,8 +147,10 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Ea
     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'
+    write(*,*)
   end do
   write(*,'(A60)')           '-------------------------------------------------'
+  write(*,*)
 
   do iEns=2,nEns
     write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy  1 ->',iEns,': ',(Omaux(iEns)+OmxcDD(iEns))*HaToeV,' eV'
@@ -157,6 +159,7 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Ea
     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'
+    write(*,*)
   end do
   write(*,'(A60)')           '-------------------------------------------------'
   write(*,*)
@@ -179,6 +182,7 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Ea
     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'
+    write(*,*)
   end do
   write(*,'(A60)')           '-------------------------------------------------'
 
@@ -192,6 +196,7 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Ea
     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'
+    write(*,*)
   end do
   write(*,'(A60)')           '-------------------------------------------------'
   write(*,*)

src/eDFT/restricted_auxiliary_energy.f90

@@ -27,7 +27,7 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
 
   Eaux(iEns) = 2d0*sum(eps(1:nO))
 
-! Doubly-excited state density matrix
+! Singly-excited state density matrix
 
   iEns = 2 
   
@@ -37,6 +37,18 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
     Eaux(iEns) = 0d0
   end if
 
+  Eaux(iEns) = Eaux(iEns) + eps(nO) + eps(nO+1)
+
+! Doubly-excited state density matrix
+
+  iEns = 3 
+  
+  if(nO > 1) then 
+    Eaux(iEns) = 2d0*sum(eps(1:nO-1))
+  else
+    Eaux(iEns) = 0d0
+  end if
+
   Eaux(iEns) = Eaux(iEns) + 2d0*eps(nO+1)
 
 end subroutine restricted_auxiliary_energy

src/eDFT/restricted_density_matrix.f90

@@ -37,6 +37,19 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
     P(:,:,iEns) = 0d0
   end if
 
+  P(:,:,iEns) = P(:,:,iEns) + 1d0*matmul(c(:,nO  :nO  ),transpose(c(:,nO  :nO  ))) &
+                            + 1d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
+
+! Doubly-excited state density matrix
+
+  iEns = 3 
+  
+  if(nO > 1) then 
+    P(:,:,iEns) = 2d0*matmul(c(:,1:nO-1),transpose(c(:,1:nO-1))) 
+  else
+    P(:,:,iEns) = 0d0
+  end if
+
   P(:,:,iEns) = P(:,:,iEns) + 2d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
 
 end subroutine restricted_density_matrix