diff --git a/input/basis b/input/basis
index cd6cb22..fb05e68 100644
--- a/input/basis
+++ b/input/basis
@@ -1,25 +1,18 @@
-1  10
-S   4
-  1    528.5000000              0.0009400        
-  2     79.3100000              0.0072140        
-  3     18.0500000              0.0359750        
-  4      5.0850000              0.1277820        
+1   3
+S   3
+  1     13.0100000              0.0196850        
+  2      1.9620000              0.1379770        
+  3      0.4446000              0.4781480        
 S   1
-  1      1.6090000              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.5363000              1.0000000        
-S   1
-  1      0.1833000              1.0000000        
+  1      0.1220000              1.0000000        
 P   1
-  1      5.9940000              1.0000000        
-P   1
-  1      1.7450000              1.0000000        
-P   1
-  1      0.5600000              1.0000000        
-D   1
-  1      4.2990000              1.0000000        
-D   1
-  1      1.2230000              1.0000000        
-F   1
-  1      2.6800000              1.0000000
-
+  1      0.7270000              1.0000000
diff --git a/input/dft b/input/dft
index deefae4..985b124 100644
--- a/input/dft
+++ b/input/dft
@@ -19,9 +19,16 @@
 # Number of states in ensemble (nEns)
   3
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
-  0.00 0.00 
+  0.00 1.00 
 # Parameters for CC weight-dependent exchange functional
-0.420431 0.069097 -0.295049
-0.135075 -0.00770826 -0.028057
+-0.00201219 -0.00371002 -0.00212719
+-0.00117715 0.00188738 -0.000414761
+# occupation numbers of orbitals nO and nO+1
+ 1.00 0.00 
+ 1.00 0.00
+ 0.00 0.00
+ 1.00 0.00
+ 1.00 0.00
+ 1.00 1.00
 # GOK-DFT: maxSCF thresh   DIIS n_diis guess_type ortho_type
-           500    0.00001   T     5      1          1
+           1000    0.00001   T     5      1          1
diff --git a/input/molecule b/input/molecule
index c78e87e..81c624a 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,4 +1,5 @@
 # nAt nEla nElb nCore nRyd
-  1   1     1   0     0
+ 2   1     1   0     0
 # Znuc   x            y           z
-  He      0.0          0.0         0.0
+  H     0.           0.         0.
+  H     0.           0.         1.4
diff --git a/input/molecule.xyz b/input/molecule.xyz
index 797b5fc..6edc99d 100644
--- a/input/molecule.xyz
+++ b/input/molecule.xyz
@@ -1,3 +1,4 @@
-  1
+  2
 
- He     0.0000000000    0.0000000000    0.0000000000
+ H      0.0000000000    0.0000000000    0.0000000000
+ H      0.0000000000    0.0000000000    0.7408481486
diff --git a/input/weight b/input/weight
index cd6cb22..fb05e68 100644
--- a/input/weight
+++ b/input/weight
@@ -1,25 +1,18 @@
-1  10
-S   4
-  1    528.5000000              0.0009400        
-  2     79.3100000              0.0072140        
-  3     18.0500000              0.0359750        
-  4      5.0850000              0.1277820        
+1   3
+S   3
+  1     13.0100000              0.0196850        
+  2      1.9620000              0.1379770        
+  3      0.4446000              0.4781480        
 S   1
-  1      1.6090000              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.5363000              1.0000000        
-S   1
-  1      0.1833000              1.0000000        
+  1      0.1220000              1.0000000        
 P   1
-  1      5.9940000              1.0000000        
-P   1
-  1      1.7450000              1.0000000        
-P   1
-  1      0.5600000              1.0000000        
-D   1
-  1      4.2990000              1.0000000        
-D   1
-  1      1.2230000              1.0000000        
-F   1
-  1      2.6800000              1.0000000
-
+  1      0.7270000              1.0000000
diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90
index 3f657a0..d1c32f7 100644
--- a/src/eDFT/GOK_RKS.f90
+++ b/src/eDFT/GOK_RKS.f90
@@ -1,5 +1,5 @@
 subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & 
-                   max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,c)
+                   max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum)
 
 ! Perform restricted Kohn-Sham calculation for ensembles
 
@@ -37,6 +37,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_
   double precision,intent(in)   :: ENuc
 
   double precision,intent(inout):: c(nBas,nBas)
+  double precision,intent(inout),dimension(2,2,3):: occnum  
 
 ! Local variables
 
@@ -189,7 +190,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_
 !   Compute density matrix 
 !------------------------------------------------------------------------
 
-    call restricted_density_matrix(nBas,nEns,nO,c(:,:),P(:,:,:))
+    call restricted_density_matrix(nBas,nEns,nO,c(:,:),P(:,:,:),occnum)
 
 !   Weight-dependent density matrix
     
@@ -344,6 +345,6 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_
 
   call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:), &
                                     nBas,nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc,eps(:),Pw(:,:),rhow(:),drhow(:,:), &
-                                    J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,E(:),Om(:))
+                                    J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,E(:),Om(:),occnum)
 
 end subroutine GOK_RKS
diff --git a/src/eDFT/GOK_UKS.f90 b/src/eDFT/GOK_UKS.f90
index 0ee5c86..1fd9ee4 100644
--- a/src/eDFT/GOK_UKS.f90
+++ b/src/eDFT/GOK_UKS.f90
@@ -1,5 +1,5 @@
 subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, &
-                   nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
+                   nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum)
 
 ! Perform unrestricted Kohn-Sham calculation for ensembles
 
@@ -30,10 +30,12 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
   double precision,intent(in)   :: X(nBas,nBas) 
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ENuc
+  double precision,intent(in),dimension(2,2,3)   :: occnum
 
 ! Local variables
 
   integer                       :: xc_rung
+  logical                       :: LDA_centered = .false.  
   integer                       :: nSCF,nBasSq
   integer                       :: n_diis
   double precision              :: conv
@@ -127,7 +129,9 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
   if(guess_type == 1) then
 
     do ispin=1,nspin
-      F(:,:,ispin) = Hc(:,:)
+      cp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
+      call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin))
+      c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
     end do
 
   else if(guess_type == 2) then
@@ -135,6 +139,10 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
     do ispin=1,nspin
       call random_number(F(:,:,ispin))
     end do
+  else
+
+    print*,'Wrong guess option'
+  stop
 
   end if
 
@@ -171,41 +179,22 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
 
     nSCF = nSCF + 1
 
-!  Transform Fock matrix in orthogonal basis
-
-    do ispin=1,nspin
-      Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
-    end do
-
-!  Diagonalize Fock matrix to get eigenvectors and eigenvalues
-
-    cp(:,:,:) = Fp(:,:,:)
-    do ispin=1,nspin
-      call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin))
-    end do
-    
-!   Back-transform eigenvectors in non-orthogonal basis
-
-    do ispin=1,nspin
-      c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
-    end do
-
-!------------------------------------------------------------------------
-!   Compute density matrix 
+!------------------------------------------------------------------------ 
+!   Compute density matrix                                               
 !------------------------------------------------------------------------
 
-    call unrestricted_density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:))
- 
-!   Weight-dependent density matrix
-    
-    Pw(:,:,:) = 0d0
-    do iEns=1,nEns
+    call unrestricted_density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:),occnum)  
+
+!   Weight-dependent density matrix 
+
+    Pw(:,:,:) = 0d0 
+    do iEns=1,nEns  
       Pw(:,:,:) = Pw(:,:,:) + wEns(iEns)*P(:,:,:,iEns)
-    end do
+    end do 
 
-!------------------------------------------------------------------------
-!   Compute one-electron density and its gradient if necessary
-!------------------------------------------------------------------------
+!------------------------------------------------------------------------ 
+!   Compute one-electron density and its gradient if necessary 
+!------------------------------------------------------------------------  
 
     do ispin=1,nspin
       do iEns=1,nEns
@@ -248,9 +237,9 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
 !   Compute exchange potential
 
     do ispin=1,nspin
-      call unrestricted_exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,nBas,Pw(:,:,ispin),   &
-                                           ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin),Fx(:,:,ispin), & 
-                                           FxHF(:,:,ispin))
+      call unrestricted_exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,&
+                                           Pw(:,:,ispin),ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin), &
+                                           Fx(:,:,ispin),FxHF(:,:,ispin))
     end do
 
 !   Compute correlation potential
@@ -282,6 +271,25 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
 
     if(minval(rcond(:)) < 1d-15) n_diis = 0
 
+!  Transform Fock matrix in orthogonal basis 
+
+    do ispin=1,nspin  
+      Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
+    end do   
+
+!  Diagonalize Fock matrix to get eigenvectors and eigenvalues  
+
+    cp(:,:,:) = Fp(:,:,:)   
+    do ispin=1,nspin  
+      call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin)) 
+    end do 
+
+!   Back-transform eigenvectors in non-orthogonal basis 
+
+    do ispin=1,nspin
+      c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
+    end do  
+
 !------------------------------------------------------------------------
 !   Compute KS energy
 !------------------------------------------------------------------------
@@ -301,19 +309,20 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
 !  Coulomb energy
 
     EJ(1) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1)))
-    EJ(2) = trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2)))
+    EJ(2) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) &
+          + 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) 
     EJ(3) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
 
 !   Exchange energy
 
     do ispin=1,nspin
-      call unrestricted_exchange_energy(x_rung,x_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, &
-                                        Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin))
+      call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &  
+                                       Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin))     
     end do
 
 !   Correlation energy
 
-    call unrestricted_correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:),Ec)
+    call unrestricted_correlation_energy(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec)
 
 !   Total energy
 
@@ -322,7 +331,7 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
 !   Check the grid accuracy by computing the number of electrons 
 
     do ispin=1,nspin
-      nEl(ispin) = electron_number(nGrid,weight(:),rhow(:,ispin))
+      nEl(ispin) = electron_number(nGrid,weight,rhow(:,ispin))
     end do
 
 !   Dump results
@@ -352,15 +361,14 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
 
 ! Compute final KS energy
 
-  call print_UKS(nBas,nO(:),eps(:,:),c(:,:,:),ENuc,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Ew)
+  call print_UKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew)
 
 !------------------------------------------------------------------------
 ! Compute individual energies from ensemble energy
 !------------------------------------------------------------------------
 
-  call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,     &
-                                      AO(:,:),dAO(:,:,:),nO(:),nV(:),T(:,:),V(:,:),ERI(:,:,:,:),ENuc,  & 
-                                      Pw(:,:,:),rhow(:,:),drhow(:,:,:),J(:,:,:),Fx(:,:,:),FxHF(:,:,:), & 
-                                      Fc(:,:,:),P(:,:,:,:),rho(:,:,:),drho(:,:,:,:),E(:),Om(:))
+
+ call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
+                                   AO,dAO,nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum) 
 
 end subroutine GOK_UKS
diff --git a/src/eDFT/LIM_RKS.f90 b/src/eDFT/LIM_RKS.f90
index a3ce34b..85b8f70 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,LDA_centered,nEns,aCC_w1,aCC_w2,nGrid,weight, &
-                   maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c)
+                   maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c,occnum)
 
 ! Perform restricted Kohn-Sham calculation for ensembles
 
@@ -30,6 +30,8 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
   double precision,intent(in)   :: X(nBas,nBas) 
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ENuc
+  double precision,intent(in),dimension(2,2,3)   :: occnum
+
 
   double precision,intent(out)  :: c(nBas,nBas) 
 
@@ -73,7 +75,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
   write(*,*)
 
   call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
-               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c)
+               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c,occnum)
 
 !------------------------------------------------------------------------
 ! Equiensemble calculation
@@ -94,7 +96,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
   write(*,*)
 
   call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
-               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c)
+               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c,occnum)
 
 !------------------------------------------------------------------------
 ! Equiensemble calculation
@@ -115,7 +117,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
   write(*,*)
 
 ! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
-!              max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c)
+!              max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c,occnum)
 
 
 !------------------------------------------------------------------------
diff --git a/src/eDFT/MOM_RKS.f90 b/src/eDFT/MOM_RKS.f90
index 01af545..cd7e00e 100644
--- a/src/eDFT/MOM_RKS.f90
+++ b/src/eDFT/MOM_RKS.f90
@@ -1,5 +1,5 @@
 subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
-                   aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c)
+                   aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c,occnum)
 
 ! Perform restricted Kohn-Sham calculation for ensembles
 
@@ -30,6 +30,8 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
   double precision,intent(in)   :: X(nBas,nBas) 
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ENuc
+  double precision,intent(in),dimension(2,2,3)   :: occnum
+
 
   double precision,intent(out)  :: c(nBas,nBas) 
 
@@ -73,7 +75,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
   write(*,*)
 
   call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
-               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c)
+               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c,occnum)
 
 !------------------------------------------------------------------------
 ! Equiensemble calculation
@@ -94,7 +96,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
   write(*,*)
 
 ! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
-!              max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c)
+!              max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c,occnum)
 
 !------------------------------------------------------------------------
 ! Equiensemble calculation
@@ -115,7 +117,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
   write(*,*)
 
   call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, &
-               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c)
+               max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c,occnum)
 
 !------------------------------------------------------------------------
 ! MOM excitation energies
diff --git a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90
index 17a92ae..caf6122 100644
--- a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90
+++ b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90
@@ -82,16 +82,21 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
 
 ! Double weight-dependency
 
-  dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) &
-         * (1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2))
+!  dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) &
+!         * (1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2))
 
-  dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2))                            &
-         * (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))  
-
-! Single weight-dependency
+!  dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2))                            &
+!         * (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))  
 
+! left single-weight-dependency
 !  dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0)))
-!  dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
+!  dCxdw2 = 0.d0
+
+! right single-weight-dependency
+  dCxdw1 = 0.d0
+  dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
+
+
 
   dCxdw1 = alpha*dCxdw1
   dCxdw2 = alpha*dCxdw2
diff --git a/src/eDFT/UCC_lda_exchange_energy.f90 b/src/eDFT/UCC_lda_exchange_energy.f90
index ad0645c..3618298 100644
--- a/src/eDFT/UCC_lda_exchange_energy.f90
+++ b/src/eDFT/UCC_lda_exchange_energy.f90
@@ -77,7 +77,14 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
   w2 = wEns(3)
   Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
 
-  Cx = alpha*Fx2*Fx1
+! for two-weights ensemble
+!  Cx = alpha*Fx2*Fx1
+
+! for left ensemble
+!   Cx = alpha*Fx1
+
+! for right ensemble
+   Cx = alpha*Fx2
 
 ! Compute GIC-LDA exchange energy
 
diff --git a/src/eDFT/UCC_lda_exchange_individual_energy.f90 b/src/eDFT/UCC_lda_exchange_individual_energy.f90
index 76560f6..25fe2cd 100644
--- a/src/eDFT/UCC_lda_exchange_individual_energy.f90
+++ b/src/eDFT/UCC_lda_exchange_individual_energy.f90
@@ -76,7 +76,14 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
   w2 = wEns(3)
   Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
 
-  Cx = alpha*Fx1*Fx2
+! for two-weight ensembles
+!  Cx = alpha*Fx1*Fx2
+
+! for left ensembles
+!  Cx = alpha*Fx1
+
+! for right ensembles
+  Cx = alpha*Fx2  
 
 ! Compute LDA exchange matrix in the AO basis
 
diff --git a/src/eDFT/UCC_lda_exchange_potential.f90 b/src/eDFT/UCC_lda_exchange_potential.f90
index 7c853e8..2dce16a 100644
--- a/src/eDFT/UCC_lda_exchange_potential.f90
+++ b/src/eDFT/UCC_lda_exchange_potential.f90
@@ -79,7 +79,14 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
   w2 = wEns(3)
   Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
 
-  Cx = alpha*Fx2*Fx1
+! for two-weight ensembles
+!  Cx = alpha*Fx2*Fx1
+
+! for left ensemble
+!   Cx = alpha*Fx1
+
+! for right ensemble
+  Cx = alpha*Fx2
 
 ! Compute LDA exchange matrix in the AO basis
 
diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90
index 46af2ae..857c733 100644
--- a/src/eDFT/eDFT.f90
+++ b/src/eDFT/eDFT.f90
@@ -7,7 +7,7 @@ program eDFT
 
   integer                       :: nNuc,nBas
   integer                       :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin)
-  double precision              :: ENuc,Ew
+  double precision              :: ENuc,Ew,ncent
 
   double precision,allocatable  :: ZNuc(:),rNuc(:,:)
 
@@ -59,6 +59,8 @@ program eDFT
   integer                       :: guess_type
   integer                       :: ortho_type
 
+  double precision,dimension(2,2,3) :: occnum  
+
 ! Hello World
 
   write(*,*)
@@ -71,7 +73,7 @@ program eDFT
 ! Read input information
 !------------------------------------------------------------------------
 
-! Read number of atoms, number of electroes of the system
+! Read number of atoms, number of electrons of the system
 ! nC   = number of core orbitals
 ! nO   = number of occupied orbitals
 ! nV   = number of virtual orbitals (see below)
@@ -82,6 +84,12 @@ program eDFT
   call read_molecule(nNuc,nEl(:),nO(:),nC(:),nR(:))
   allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
 
+  ncent = dble(nEl(1) + nEl(2))
+  print*, 'ncent=',ncent
+  print*, 'N-1/N=',(ncent-1.d0)/ncent
+  print*, 'N+1/N=',(ncent+1.d0)/ncent
+
+
 ! Read geometry
 
   call read_geometry(nNuc,ZNuc,rNuc,ENuc)
@@ -100,11 +108,12 @@ program eDFT
 ! DFT options
 !------------------------------------------------------------------------
 
+
 ! Allocate ensemble weights
 
   allocate(wEns(maxEns))
   call read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, & 
-                    maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type)
+                    maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,ncent,occnum)
 
 !------------------------------------------------------------------------
 ! Read one- and two-electron integrals
@@ -166,7 +175,7 @@ program eDFT
     call cpu_time(start_KS)
     call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, &
                  maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO(1),nV(1), &
-                 S,T,V,Hc,ERI,X,ENuc,Ew,c)
+                 S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum)
     call cpu_time(end_KS)
 
     t_KS = end_KS - start_KS
@@ -184,7 +193,7 @@ program eDFT
     call cpu_time(start_KS)
     call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:),           &
                  aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & 
-                 S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:))
+                 S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum)
     call cpu_time(end_KS)
 
     t_KS = end_KS - start_KS
@@ -202,7 +211,7 @@ program eDFT
     call cpu_time(start_KS)
     call MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:),           &
                  aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & 
-                 S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:))
+                 S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum)
     call cpu_time(end_KS)
 
     t_KS = end_KS - start_KS
@@ -219,7 +228,7 @@ program eDFT
 
     call cpu_time(start_KS)
     call GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, & 
-                   nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew)
+                   nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum)
     call cpu_time(end_KS)
 
     t_KS = end_KS - start_KS
@@ -236,7 +245,7 @@ program eDFT
 
     call cpu_time(start_KS)
     call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & 
-                  nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
+                  nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum)
     call cpu_time(end_KS)
 
     t_KS = end_KS - start_KS
diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90
index 1ba8e6c..dd66abd 100644
--- a/src/eDFT/eDFT_UKS.f90
+++ b/src/eDFT/eDFT_UKS.f90
@@ -1,5 +1,5 @@
 subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh,max_diis,guess_type, &
-                   nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
+                   nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum)
 
 ! Perform unrestricted Kohn-Sham calculation for ensembles
 
@@ -30,6 +30,8 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
   double precision,intent(in)   :: X(nBas,nBas) 
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ENuc
+  double precision,intent(in),dimension(2,2,3)   :: occnum
+
 
 ! Local variables
 
@@ -184,7 +186,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
 !   Compute density matrix 
 !------------------------------------------------------------------------
 
-    call unrestricted_density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:))
+    call unrestricted_density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:),occnum)
  
 !   Weight-dependent density matrix
     
@@ -369,6 +371,6 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
 !------------------------------------------------------------------------
 
   call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
-                                      AO,dAO,nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om)
+                                      AO,dAO,nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum)
 
 end subroutine eDFT_UKS
diff --git a/src/eDFT/read_options.f90 b/src/eDFT/read_options.f90
index f133ba0..e92eb7f 100644
--- a/src/eDFT/read_options.f90
+++ b/src/eDFT/read_options.f90
@@ -1,5 +1,5 @@
 subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, & 
-                        maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type)
+                        maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,ncent,occnum)
 
 ! Read DFT options
 
@@ -9,7 +9,7 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
 
 ! Local variables
 
-  integer                       :: I
+  integer                       :: I,J
 
 ! Output variables
 
@@ -21,6 +21,7 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
   double precision,intent(out)  :: wEns(maxEns)
   double precision,intent(out)  :: aCC_w1(3)
   double precision,intent(out)  :: aCC_w2(3)
+  double precision,intent(out),dimension(2,2,3)  :: occnum
 
   integer,intent(out)           :: maxSCF
   double precision,intent(out)  :: thresh
@@ -28,6 +29,7 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
   integer,intent(out)           :: max_diis
   integer,intent(out)           :: guess_type
   integer,intent(out)           :: ortho_type
+  double precision,intent(in)   :: ncent
 
 ! Local variables
 
@@ -92,10 +94,18 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
   write(*,*)'----------------------------------------------------------'
   write(*,*) 
   
-! Read ensemble weights
+! Read ensemble weights for unphysical (integer number of electrons) left or right N-centered ensembles:
+! (alpha,0) or (0,alpha) in input
+!  read(1,*)
+!  read(1,*) (wEns(I),I=2,nEns)
+!  wEns(2) = (ncent/(ncent-1.d0))*wEns(2)
+!  wEns(3) = (ncent/(ncent+1.d0))*wEns(3)
+!  wEns(1) = 1d0 - ((ncent-1.d0)/ncent)*wEns(2) - ((ncent+1.d0)/ncent)*wEns(3) ! for N-centered
+
+! Read ensemble weights for real physical (fractional number of electrons) ensemble (w1,w2)
   read(1,*)
   read(1,*) (wEns(I),I=2,nEns)
-  wEns(1) = 1d0 - sum(wEns)
+  wEns(1) = 1d0 - wEns(2) - wEns(3) 
 
   write(*,*)'----------------------------------------------------------'
   write(*,*)' Ensemble weights '
@@ -119,6 +129,14 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
   write(*,*)'----------------------------------------------------------'
   call matout(3,1,aCC_w2)
   write(*,*) 
+ 
+!  allocate(occnum(2,2,nEns))
+! Read occupation numbers for orbitals nO and nO+1
+  read(1,*)
+  do J=1,3
+    read(1,*) (occnum(1,I,J),I=1,2)
+    read(1,*) (occnum(2,I,J),I=1,2)
+  end do
 
 ! Read KS options
 
diff --git a/src/eDFT/restricted_auxiliary_energy.f90 b/src/eDFT/restricted_auxiliary_energy.f90
index 16490f5..8c294e0 100644
--- a/src/eDFT/restricted_auxiliary_energy.f90
+++ b/src/eDFT/restricted_auxiliary_energy.f90
@@ -1,4 +1,4 @@
-subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
+subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
 
 ! Compute the auxiliary KS energies 
 
@@ -12,6 +12,8 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
   integer,intent(in)            :: nEns
   integer,intent(in)            :: nO
   double precision,intent(in)   :: eps(nBas)
+  double precision,intent(in),dimension(2,2,3)   :: occnum
+
 
 ! Local variables
 
@@ -36,8 +38,8 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
   else
     Eaux(iEns) = 0d0
   end if
+  Eaux(iEns) = Eaux(iEns) + sum(occnum(:,1,iEns))*eps(nO) + sum(occnum(:,2,iEns))*eps(nO+1)
 
-  Eaux(iEns) = Eaux(iEns) + eps(nO) + eps(nO+2)
 
 ! Doubly-excited state density matrix
 
@@ -48,7 +50,7 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
   else
     Eaux(iEns) = 0d0
   end if
+  Eaux(iEns) = Eaux(iEns) + sum(occnum(:,1,iEns))*eps(nO) + sum(occnum(:,2,iEns))*eps(nO+1)
 
-  Eaux(iEns) = Eaux(iEns) + 2d0*eps(nO+1)
 
 end subroutine restricted_auxiliary_energy
diff --git a/src/eDFT/restricted_density_matrix.f90 b/src/eDFT/restricted_density_matrix.f90
index 93ac4b4..661113e 100644
--- a/src/eDFT/restricted_density_matrix.f90
+++ b/src/eDFT/restricted_density_matrix.f90
@@ -1,4 +1,4 @@
-subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
+subroutine restricted_density_matrix(nBas,nEns,nO,c,P,occnum)
 
 ! Calculate density matrices
 
@@ -12,6 +12,8 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
   integer,intent(in)            :: nEns
   integer,intent(in)            :: nO
   double precision,intent(in)   :: c(nBas,nBas)
+  double precision,intent(in),dimension(2,2,3)   :: occnum
+
 
 ! Local variables
 
@@ -25,7 +27,7 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
 
   iEns = 1
 
-  P(:,:,iEns) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
+  P(:,:,iEns) = 2.d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
 
 ! Doubly-excited state density matrix
 
@@ -36,9 +38,8 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
   else
     P(:,:,iEns) = 0d0
   end if
-
-  P(:,:,iEns) = P(:,:,iEns) + 1d0*matmul(c(:,nO  :nO  ),transpose(c(:,nO  :nO  ))) &
-                            + 1d0*matmul(c(:,nO+2:nO+2),transpose(c(:,nO+2:nO+2)))
+  P(:,:,iEns) = P(:,:,iEns) + sum(occnum(:,1,iEns))*matmul(c(:,nO:nO),transpose(c(:,nO:nO))) &
+                            + sum(occnum(:,2,iEns))*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
 
 ! Doubly-excited state density matrix
 
@@ -49,7 +50,7 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
   else
     P(:,:,iEns) = 0d0
   end if
-
-  P(:,:,iEns) = P(:,:,iEns) + 2d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
+  P(:,:,iEns) = P(:,:,iEns) + sum(occnum(:,1,iEns))*matmul(c(:,nO:nO),transpose(c(:,nO:nO))) &
+              + sum(occnum(:,2,iEns))*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
 
 end subroutine restricted_density_matrix
diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90
index 67968f0..466e300 100644
--- a/src/eDFT/restricted_individual_energy.f90
+++ b/src/eDFT/restricted_individual_energy.f90
@@ -1,5 +1,5 @@
 subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, & 
-                                        nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,E,Om)
+                                        nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,E,Om,occnum)
 
 ! Compute individual energies as well as excitation energies
 
@@ -38,6 +38,8 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n
   double precision,intent(in)   :: J(nBas,nBas)
 
   double precision              :: Ew
+  double precision,intent(in),dimension(2,2,3)   :: occnum
+
 
 ! Local variables
 
@@ -108,7 +110,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n
 ! Compute auxiliary energies
 !------------------------------------------------------------------------
 
-  call restricted_auxiliary_energy(nBas,nEns,nO,eps(:),Eaux(:))
+  call restricted_auxiliary_energy(nBas,nEns,nO,eps(:),Eaux(:),occnum)
 
 !------------------------------------------------------------------------
 ! Compute derivative discontinuities
diff --git a/src/eDFT/unrestricted_auxiliary_energy.f90 b/src/eDFT/unrestricted_auxiliary_energy.f90
index 192f6dd..be14e03 100644
--- a/src/eDFT/unrestricted_auxiliary_energy.f90
+++ b/src/eDFT/unrestricted_auxiliary_energy.f90
@@ -1,4 +1,4 @@
-subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
+subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
 
 ! Compute the auxiliary KS energies 
 
@@ -12,6 +12,8 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
   integer,intent(in)            :: nEns
   integer,intent(in)            :: nO(nspin)
   double precision,intent(in)   :: eps(nBas,nspin)
+  double precision,intent(in),dimension(2,2,3)   :: occnum
+
 
 ! Local variables
 
@@ -21,30 +23,74 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
 
   double precision,intent(out)  :: Eaux(nspin,nEns)
 
-! N-electron ground state
-
+!----------------------------------------------------------
+!-------------- GOK-UKS ----------------------------------
+!----------------------------------------------------------
   iEns = 1
   do ispin=1,nspin
     Eaux(ispin,iEns) = sum(eps(1:nO(ispin),ispin))
   end do
 
-! (N-1)-electron ground state
-
-  iEns = 2 
- 
-  Eaux(2,iEns) = sum(eps(1:nO(2),2))
-
-  if(nO(1) > 1) then 
+  iEns = 2
+  if(nO(2) > 1) then
+    Eaux(2,iEns) = sum(eps(1:nO(2)-1,2))
+  else
+    Eaux(2,iEns) = 0d0
+  end if
+  Eaux(2,iEns) = sum(eps(1:nO(2)-1,2)) + occnum(2,1,iEns)*eps(nO(2),2) &
+               + occnum(2,2,iEns)*eps(nO(2)+1,2)
+  if(nO(1) > 1) then
     Eaux(1,iEns) = sum(eps(1:nO(1)-1,1))
   else
     Eaux(1,iEns) = 0d0
   end if
+  Eaux(1,iEns) = Eaux(1,iEns) + occnum(1,1,iEns)*eps(nO(1),1) &
+               + occnum(1,2,iEns)*eps(nO(1)+1,1)
+
+  iEns = 3
+  if(nO(1) > 1) then
+    Eaux(1,iEns) = sum(eps(1:nO(1)-1,1))
+  else
+    Eaux(1,iEns) = 0d0
+  end if
+  Eaux(1,iEns) = Eaux(1,iEns) + occnum(1,1,iEns)*eps(nO(1),1) &
+               + occnum(1,2,iEns)*eps(nO(1)+1,1)
+
+  if(nO(2) > 1) then
+    Eaux(2,iEns) = sum(eps(1:nO(2)-1,2))
+  else
+    Eaux(2,iEns) = 0d0
+  end if
+  Eaux(2,iEns) = Eaux(2,iEns) + occnum(2,1,iEns)*eps(nO(2),2) &
+               + occnum(2,2,iEns)* eps(nO(2)+1,2)
+
+!----------------------------------------------------------
+!-------------- eDFT-UKS ----------------------------------
+!----------------------------------------------------------
+! N-electron ground state
+
+!  iEns = 1
+!  do ispin=1,nspin
+!    Eaux(ispin,iEns) = sum(eps(1:nO(ispin),ispin))
+!  end do
+
+! (N-1)-electron ground state
+
+!  iEns = 2 
+ 
+!  Eaux(2,iEns) = sum(eps(1:nO(2),2))
+
+!  if(nO(1) > 1) then 
+!    Eaux(1,iEns) = sum(eps(1:nO(1)-1,1))
+!  else
+!    Eaux(1,iEns) = 0d0
+!  end if
 
 ! (N+1)-electron ground state
 
-  iEns = 3 
+!  iEns = 3 
  
-  Eaux(2,iEns) = sum(eps(1:nO(2)+1,2))
-  Eaux(1,iEns) = sum(eps(1:nO(1),1)) 
+!  Eaux(2,iEns) = sum(eps(1:nO(2)+1,2))
+!  Eaux(1,iEns) = sum(eps(1:nO(1),1)) 
 
 end subroutine unrestricted_auxiliary_energy
diff --git a/src/eDFT/unrestricted_density_matrix.f90 b/src/eDFT/unrestricted_density_matrix.f90
index 372a063..36de196 100644
--- a/src/eDFT/unrestricted_density_matrix.f90
+++ b/src/eDFT/unrestricted_density_matrix.f90
@@ -1,4 +1,4 @@
-subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P)
+subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P,occnum)
 
 ! Calculate density matrices
 
@@ -12,6 +12,8 @@ subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P)
   integer,intent(in)            :: nEns
   integer,intent(in)            :: nO(nspin)
   double precision,intent(in)   :: c(nBas,nBas,nspin)
+  double precision,intent(in),dimension(2,2,3)   :: occnum
+
 
 ! Local variables
 
@@ -22,26 +24,75 @@ subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P)
 
   double precision,intent(out)  :: P(nBas,nBas,nspin,nEns)
 
-! N-electron ground state
-
+!-------------------------------------------------------
+!------------------------ GOK-UKS ----------------------
+!-------------------------------------------------------
   iEns = 1
   do ispin=1,nspin
-    P(:,:,ispin,iEns) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
+    P(:,:,ispin,iEns) =  matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
   end do
 
+  iEns = 2
+
+  if(nO(1) > 1) then 
+    P(:,:,1,iEns) = matmul(c(:,1:nO(1)-1,1),transpose(c(:,1:nO(1)-1,1)))
+  else
+    P(:,:,1,iEns)=0.d0 
+  end if
+  P(:,:,1,iEns) = P(:,:,1,iEns) + occnum(1,1,iEns)* matmul(c(:,nO(1):nO(1),1),transpose(c(:,nO(1):nO(1),1))) &
+                + occnum(1,2,iEns)* matmul(c(:,nO(1)+1:nO(1)+1,1),transpose(c(:,nO(1)+1:nO(1)+1,1)))
+ if(nO(2) > 1) then
+    P(:,:,2,iEns) = matmul(c(:,1:nO(2)-1,2),transpose(c(:,1:nO(2)-1,2)))
+  else
+    P(:,:,2,iEns)=0.d0
+  end if
+
+  P(:,:,2,iEns) = P(:,:,2,iEns) + occnum(2,1,iEns)* matmul(c(:,1:nO(2),2),transpose(c(:,1:nO(2),2))) &
+                + occnum(2,2,iEns)*matmul(c(:,1:nO(2)+1,2),transpose(c(:,1:nO(2)+1,2)))
+
+ 
+  iEns = 3  
+
+  if(nO(1) > 1) then
+    P(:,:,1,iEns) = matmul(c(:,1:nO(1)-1,1),transpose(c(:,1:nO(1)-1,1)))
+  else
+    P(:,:,1,iEns)=0.d0
+  end if
+  P(:,:,1,iEns) = P(:,:,1,iEns) + occnum(1,1,iEns)* matmul(c(:,nO(1):nO(1),1),transpose(c(:,nO(1):nO(1),1))) &
+                + occnum(1,2,iEns)* matmul(c(:,nO(1)+1:nO(1)+1,1),transpose(c(:,nO(1)+1:nO(1)+1,1)))
+
+  if(nO(2) > 1) then
+    P(:,:,2,iEns) = matmul(c(:,1:nO(2)-1,2),transpose(c(:,1:nO(2)-1,2)))
+  else
+    P(:,:,2,iEns)=0.d0
+  end if
+  P(:,:,2,iEns) = P(:,:,2,iEns) + occnum(2,1,iEns)* matmul(c(:,nO(2):nO(2),2),transpose(c(:,nO(2):nO(2),2))) &
+                + occnum(2,2,iEns)*matmul(c(:,nO(2)+1:nO(2)+1,2),transpose(c(:,nO(2)+1:nO(2)+1,2)))
+
+!-------------------------------------------------------------------
+!--------------- For eDFT_UKS (fundamental gap)---------------------
+!-------------------------------------------------------------------
+
+! N-electron ground state
+
+!  iEns = 1
+!  do ispin=1,nspin
+!    P(:,:,ispin,iEns) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
+!  end do
+
 ! (N-1)-electron state: remove spin-up electrons
 
-  iEns = 2
-  P(:,:,2,iEns) = matmul(c(:,1:nO(2),2),transpose(c(:,1:nO(2),2)))  
-  if (nO(1) > 1) then
-    P(:,:,1,iEns) = matmul(c(:,1:nO(1)-1,1),transpose(c(:,1:nO(1)-1,1)))
-  else
-    P(:,:,1,iEns) = 0.d0
-  end if
+!  iEns = 2
+!  P(:,:,2,iEns) = matmul(c(:,1:nO(2),2),transpose(c(:,1:nO(2),2)))  
+!  if (nO(1) > 1) then
+!    P(:,:,1,iEns) = matmul(c(:,1:nO(1)-1,1),transpose(c(:,1:nO(1)-1,1)))
+!  else
+!    P(:,:,1,iEns) = 0.d0
+!  end if
 ! (N+1)-electron state: remove spin-down electrons
 
-  iEns = 3 
-  P(:,:,2,iEns) = matmul(c(:,1:nO(2)+1,2),transpose(c(:,1:nO(2)+1,2)))  
-  P(:,:,1,iEns) = matmul(c(:,1:nO(1),1),transpose(c(:,1:nO(1),1)))
+!  iEns = 3 
+!  P(:,:,2,iEns) = matmul(c(:,1:nO(2)+1,2),transpose(c(:,1:nO(2)+1,2)))  
+!  P(:,:,1,iEns) = matmul(c(:,1:nO(1),1),transpose(c(:,1:nO(1),1)))
 
 end subroutine unrestricted_density_matrix
diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90
index 6655569..1e6d58d 100644
--- a/src/eDFT/unrestricted_individual_energy.f90
+++ b/src/eDFT/unrestricted_individual_energy.f90
@@ -1,5 +1,5 @@
 subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,dAO, &
-                                          nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om)
+                                          nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum)
 
 ! Compute unrestricted individual energies as well as excitation energies
 
@@ -41,6 +41,8 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
   double precision,intent(in)   :: FxHF(nBas,nBas,nspin)
   double precision,intent(in)   :: Fc(nBas,nBas,nspin)
   double precision              :: Ew
+  double precision,intent(in),dimension(2,2,3)   :: occnum
+
 
 ! Local variables
 
@@ -182,7 +184,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
 ! Compute auxiliary energies
 !------------------------------------------------------------------------
  
-  call unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
+  call unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
 
 !------------------------------------------------------------------------
 ! Compute derivative discontinuities