diff --git a/examples/basis.C2-cc-pvdz b/examples/basis.C2-cc-pvdz
index 16d91fa..f44df8f 100644
--- a/examples/basis.C2-cc-pvdz
+++ b/examples/basis.C2-cc-pvdz
@@ -1,64 +1,58 @@
-1  6
-S   9
-1         6.665000E+03           6.920000E-04
-2         1.000000E+03           5.329000E-03
-3         2.280000E+02           2.707700E-02
-4         6.471000E+01           1.017180E-01
-5         2.106000E+01           2.747400E-01
-6         7.495000E+00           4.485640E-01
-7         2.797000E+00           2.850740E-01
-8         5.215000E-01           1.520400E-02
-9         1.596000E-01          -3.191000E-03
-S   9
-1         6.665000E+03          -1.460000E-04
-2         1.000000E+03          -1.154000E-03
-3         2.280000E+02          -5.725000E-03
-4         6.471000E+01          -2.331200E-02
-5         2.106000E+01          -6.395500E-02
-6         7.495000E+00          -1.499810E-01
-7         2.797000E+00          -1.272620E-01
-8         5.215000E-01           5.445290E-01
-9         1.596000E-01           5.804960E-01
+1   6
+S   8
+  1   6665.0000000              0.0006920
+  2    1000.000000              0.0053290
+  3    228.0000000              0.0270770
+  4     64.7100000              0.1017180
+  5      21.060000              0.2747400
+  6      7.4950000              0.4485640    
+  7      2.7970000              0.2850740    
+  8      0.5215000              0.0152040    
+S   8
+  1   6665.0000000             -0.0001460
+  2    1000.000000             -0.0011540    
+  3    228.0000000             -0.0057250    
+  4     64.7100000             -0.0233120    
+  5      21.060000             -0.0639550    
+  6      7.4950000             -0.1499810    
+  7      2.7970000             -0.1272620    
+  8      0.5215000              0.5445290  
 S   1
-1         1.596000E-01           1.000000E+00
-P   4
-1         9.439000E+00           3.810900E-02
-2         2.002000E+00           2.094800E-01
-3         5.456000E-01           5.085570E-01
-4         1.517000E-01           4.688420E-01
+  1      0.1596000              1.0000000    
+P   3
+  1      9.4390000              0.0381090    
+  2      2.0020000              0.2094800
+  3      0.5456000              0.5085570
 P   1
-1         1.517000E-01           1.000000E+00
+  1      0.1517000              1.0000000
 D   1
-1         5.500000E-01           1.0000000
+  1      0.5500000              1.0000000
 2   6
-S   9
-1         6.665000E+03           6.920000E-04
-2         1.000000E+03           5.329000E-03
-3         2.280000E+02           2.707700E-02
-4         6.471000E+01           1.017180E-01
-5         2.106000E+01           2.747400E-01
-6         7.495000E+00           4.485640E-01
-7         2.797000E+00           2.850740E-01
-8         5.215000E-01           1.520400E-02
-9         1.596000E-01          -3.191000E-03
-S   9
-1         6.665000E+03          -1.460000E-04
-2         1.000000E+03          -1.154000E-03
-3         2.280000E+02          -5.725000E-03
-4         6.471000E+01          -2.331200E-02
-5         2.106000E+01          -6.395500E-02
-6         7.495000E+00          -1.499810E-01
-7         2.797000E+00          -1.272620E-01
-8         5.215000E-01           5.445290E-01
-9         1.596000E-01           5.804960E-01
+S   8
+  1   6665.0000000              0.0006920
+  2    1000.000000              0.0053290
+  3    228.0000000              0.0270770
+  4     64.7100000              0.1017180
+  5      21.060000              0.2747400
+  6      7.4950000              0.4485640    
+  7      2.7970000              0.2850740    
+  8      0.5215000              0.0152040    
+S   8
+  1   6665.0000000             -0.0001460
+  2    1000.000000             -0.0011540    
+  3    228.0000000             -0.0057250    
+  4     64.7100000             -0.0233120    
+  5      21.060000             -0.0639550    
+  6      7.4950000             -0.1499810    
+  7      2.7970000             -0.1272620    
+  8      0.5215000              0.5445290    
 S   1
-1         1.596000E-01           1.000000E+00
-P   4
-1         9.439000E+00           3.810900E-02
-2         2.002000E+00           2.094800E-01
-3         5.456000E-01           5.085570E-01
-4         1.517000E-01           4.688420E-01
+  1      0.1596000              1.0000000    
+P   3
+  1      9.4390000              0.0381090    
+  2      2.0020000              0.2094800
+  3      0.5456000              0.5085570
 P   1
-1         1.517000E-01           1.000000E+00
+  1      0.1517000              1.0000000
 D   1
-1         5.500000E-01           1.0000000
+  1      0.5500000              1.0000000
diff --git a/input/basis b/input/basis
index fb05e68..a82aa28 100644
--- a/input/basis
+++ b/input/basis
@@ -1,18 +1,29 @@
-1   3
-S   3
-  1     13.0100000              0.0196850        
-  2      1.9620000              0.1379770        
-  3      0.4446000              0.4781480        
+1   6
+S   8
+  1   1469.0000000              0.0007660        
+  2    220.5000000              0.0058920        
+  3     50.2600000              0.0296710        
+  4     14.2400000              0.1091800        
+  5      4.5810000              0.2827890        
+  6      1.5800000              0.4531230        
+  7      0.5640000              0.2747740        
+  8      0.0734500              0.0097510        
+S   8
+  1   1469.0000000             -0.0001200        
+  2    220.5000000             -0.0009230        
+  3     50.2600000             -0.0046890        
+  4     14.2400000             -0.0176820        
+  5      4.5810000             -0.0489020        
+  6      1.5800000             -0.0960090        
+  7      0.5640000             -0.1363800        
+  8      0.0734500              0.5751020        
 S   1
-  1      0.1220000              1.0000000        
+  1      0.0280500              1.0000000        
+P   3
+  1      1.5340000              0.0227840        
+  2      0.2749000              0.1391070        
+  3      0.0736200              0.5003750        
 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.1220000              1.0000000        
-P   1
-  1      0.7270000              1.0000000
+  1      0.0240300              1.0000000        
+D   1
+  1      0.1239000		1.0000000
diff --git a/input/dft b/input/dft
index 8e8e1f6..a4a3113 100644
--- a/input/dft
+++ b/input/dft
@@ -18,21 +18,23 @@
   1
 # Number of states in ensemble (nEns)
   3
+# occupation numbers of orbitals
+1 1 0 0 0 0 0 0 0 0 0 0 0 0 0  
+1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+
+1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+
+1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
-  1.00 0.00 
+  0.20 0.30 
 # Ncentered ? 0 for NO 
- 0 
+ 1 
 # Parameters for CC weight-dependent exchange functional
 0.420431 0.069097 -0.295049
 0.135075 -0.00770826 -0.028057
 # choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2
 2
-# 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
            1000    0.00001   T     5      1          1
diff --git a/input/molecule b/input/molecule
index 81c624a..058d6dd 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,5 +1,4 @@
 # nAt nEla nElb nCore nRyd
- 2   1     1   0     0
+  1   2     1   0     0
 # Znuc   x            y           z
-  H     0.           0.         0.
-  H     0.           0.         1.4
+  Li     0.0          0.0         0.0
diff --git a/input/molecule.xyz b/input/molecule.xyz
index 6edc99d..c9a5a65 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
+ Li     0.0000000000    0.0000000000    0.0000000000
diff --git a/input/weight b/input/weight
index fb05e68..a82aa28 100644
--- a/input/weight
+++ b/input/weight
@@ -1,18 +1,29 @@
-1   3
-S   3
-  1     13.0100000              0.0196850        
-  2      1.9620000              0.1379770        
-  3      0.4446000              0.4781480        
+1   6
+S   8
+  1   1469.0000000              0.0007660        
+  2    220.5000000              0.0058920        
+  3     50.2600000              0.0296710        
+  4     14.2400000              0.1091800        
+  5      4.5810000              0.2827890        
+  6      1.5800000              0.4531230        
+  7      0.5640000              0.2747740        
+  8      0.0734500              0.0097510        
+S   8
+  1   1469.0000000             -0.0001200        
+  2    220.5000000             -0.0009230        
+  3     50.2600000             -0.0046890        
+  4     14.2400000             -0.0176820        
+  5      4.5810000             -0.0489020        
+  6      1.5800000             -0.0960090        
+  7      0.5640000             -0.1363800        
+  8      0.0734500              0.5751020        
 S   1
-  1      0.1220000              1.0000000        
+  1      0.0280500              1.0000000        
+P   3
+  1      1.5340000              0.0227840        
+  2      0.2749000              0.1391070        
+  3      0.0736200              0.5003750        
 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.1220000              1.0000000        
-P   1
-  1      0.7270000              1.0000000
+  1      0.0240300              1.0000000        
+D   1
+  1      0.1239000		1.0000000
diff --git a/scan_w.sh b/scan_w.sh
index fb0583a..799d780 100755
--- a/scan_w.sh
+++ b/scan_w.sh
@@ -76,11 +76,11 @@ do
   echo "2"  >> input/dft
   echo "# occupation numbers of orbitals nO and nO+1"  >> input/dft
   echo " 1.00 0.00 " >> input/dft
-  echo " 1.00 0.00 " >> input/dft
+  echo " 0.00 0.00 " >> input/dft
+  echo " 0.00 0.00 " >> input/dft
   echo " 0.00 0.00 " >> input/dft
   echo " 1.00 0.00 " >> input/dft
-  echo " 1.00 0.00 " >> input/dft
-  echo " 1.00 1.00 " >> input/dft
+  echo " 0.00 1.00 " >> input/dft
   echo "# GOK-DFT: maxSCF thresh   DIIS n_diis guess_type ortho_type" >> input/dft
   echo "           1000    0.00001   T     5      1          1" >> input/dft
   OUTPUT=${MOL}_${BASIS}_${XF}_${CF}_${w2}.out
diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90
index 7186320..7aa201d 100644
--- a/src/eDFT/GOK_RKS.f90
+++ b/src/eDFT/GOK_RKS.f90
@@ -37,7 +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(in):: occnum(nspin,2,nEns)  
+  double precision,intent(in)   :: occnum(nBas,nspin,nEns)  
   integer,intent(in)            :: Cx_choice
 
 ! Local variables
diff --git a/src/eDFT/LIM_RKS.f90 b/src/eDFT/LIM_RKS.f90
index 6893023..263fb51 100644
--- a/src/eDFT/LIM_RKS.f90
+++ b/src/eDFT/LIM_RKS.f90
@@ -31,7 +31,7 @@ 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)   :: occnum(nspin,2,nEns)
+  double precision,intent(in)   :: occnum(nBas,nspin,nEns)
   integer,intent(in)            :: Cx_choice
 
 
diff --git a/src/eDFT/MOM_RKS.f90 b/src/eDFT/MOM_RKS.f90
index 3f54884..35e2c8c 100644
--- a/src/eDFT/MOM_RKS.f90
+++ b/src/eDFT/MOM_RKS.f90
@@ -31,7 +31,7 @@ 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)   :: occnum(nspin,2,nEns)
+  double precision,intent(in)   :: occnum(nBas,nspin,nEns)
   integer,intent(in)            :: Cx_choice
 
 
diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90
index 10d6465..e9ba259 100644
--- a/src/eDFT/eDFT.f90
+++ b/src/eDFT/eDFT.f90
@@ -85,12 +85,6 @@ 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)
@@ -109,12 +103,11 @@ program eDFT
 ! DFT options
 !------------------------------------------------------------------------
 
-
 ! Allocate ensemble weights
-  allocate(occnum(nspin,2,maxEns))
-  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,doNcentered,ncent,occnum,Cx_choice)
+
+  allocate(wEns(maxEns),occnum(nBas,nspin,maxEns))
+  call read_options(nBas,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,doNcentered,occnum,Cx_choice)
 
 !------------------------------------------------------------------------
 ! Read one- and two-electron integrals
diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90
index fddf900..49f2b87 100644
--- a/src/eDFT/eDFT_UKS.f90
+++ b/src/eDFT/eDFT_UKS.f90
@@ -30,7 +30,7 @@ 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)   :: occnum(nspin,2,nEns)
+  double precision,intent(in)   :: occnum(nBas,nspin,nEns)
   integer,intent(in)            :: Cx_choice
 
 ! Local variables
@@ -186,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(:,:,:,:),occnum)
+    call unrestricted_density_matrix(nBas,nEns,c(:,:,:),P(:,:,:,:),occnum(:,:,:))
  
 !   Weight-dependent density matrix
     
@@ -337,6 +337,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
       nEl(ispin) = electron_number(nGrid,weight,rhow(:,ispin))
     end do
 
+
 !   Dump results
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & 
@@ -362,6 +363,17 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
 
   end if
 
+!!!!!
+    do iEns=1,nEns
+      print*,'occnum=',occnum(1,1,iEns),occnum(2,1,iEns),occnum(1,2,iEns),occnum(2,2,iEns)
+      print*,'nel up and down and total=', electron_number(nGrid,weight,&
+                                           rho(:,1,iEns)),electron_number(nGrid,weight,rho(:,2,iEns)),sum(nEl(:))
+
+    end do
+!!!!!
+
+
+
 ! Compute final KS energy
 
   call print_UKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew)
diff --git a/src/eDFT/read_options.f90 b/src/eDFT/read_options.f90
index 7e6b036..4eb7f62 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,doNcentered,ncent,occnum,Cx_choice)
+subroutine read_options(nBas,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,doNcentered,occnum,Cx_choice)
 
 ! Read DFT options
 
@@ -7,9 +7,16 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
 
   include 'parameters.h'
 
+! Input variables
+  integer,intent(in)           :: nBas
+
 ! Local variables
 
-  integer                       :: I,J
+  integer                       :: iBas
+  integer                       :: iEns
+  integer                       :: iParam
+  character(len=1)              :: answer
+  double precision,allocatable  :: nEl(:)
 
 ! Output variables
 
@@ -17,11 +24,12 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
   integer,intent(out)           :: x_rung,c_rung
   character(len=12),intent(out) :: x_DFA, c_DFA
   integer,intent(out)           :: SGn
-  integer,intent(out)           :: nEns, doNcentered
+  integer,intent(out)           :: nEns
+  integer,intent(out)           :: doNcentered
   double precision,intent(out)  :: wEns(maxEns)
   double precision,intent(out)  :: aCC_w1(3)
   double precision,intent(out)  :: aCC_w2(3)
-  double precision,intent(inout):: occnum(nspin,2,maxEns)
+  double precision,intent(out)  :: occnum(nBas,nspin,maxEns)
 
   integer,intent(out)           :: maxSCF
   double precision,intent(out)  :: thresh
@@ -29,13 +37,8 @@ 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
   integer,intent(out)           :: Cx_choice
 
-! Local variables
-
-  character(len=1)              :: answer
-
 ! Open file with method specification
 
   open(unit=1,file='input/dft')
@@ -95,26 +98,60 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
   write(*,*)'----------------------------------------------------------'
   write(*,*) 
   
-! 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 occupation numbers for orbitals nO and nO+1
 
+  occnum(:,:,:) = 0d0
+
+  do iEns=1,nEns
+    read(1,*)
+    read(1,*) (occnum(iBas,1,iEns),iBas=1,nBas)
+    read(1,*) (occnum(iBas,2,iEns),iBas=1,nBas)
+  end do
+
+  do iEns=1,nEns
+    write(*,*) 
+    write(*,*) '==============='
+    write(*,*) 'State n.',iEns
+    write(*,*) '==============='
+    write(*,*)
+    write(*,*) 'Spin-up   occupation numbers'
+    write(*,*) (int(occnum(iBas,1,iEns)),iBas=1,nBas)
+    write(*,*) 'Spin-down occupation numbers'
+    write(*,*) (int(occnum(iBas,2,iEns)),iBas=1,nBas)
+    write(*,*)
+  end do
 ! Read ensemble weights for real physical (fractional number of electrons) ensemble (w1,w2)
+
+  allocate(nEl(nEns))
+  nEl(:) = 0d0
+  do iEns=1,nEns
+    do iBas=1,nBas
+      nEl(iEns) = nEl(iEns) + occnum(iBas,1,iEns) + occnum(iBas,2,iEns)
+    end do
+  end do
+  print*,'nEl'
+  print*,nEl
+  
+
   read(1,*)
-  read(1,*) (wEns(I),I=2,nEns)
+  read(1,*) (wEns(iEns),iEns=2,nEns)
   read(1,*)
   read(1,*) doNcentered
+
   if (doNcentered==0) then
+
     wEns(1) = 1d0 - wEns(2) - wEns(3) 
+
   else
-  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
+
+!   wEns(1) = 1d0 -  nEl(2)/nEl(1)*wEns(2) - nEl(3)/nEl(1)*wEns(3) 
+
+    wEns(1) = 1d0 - wEns(2) - wEns(3)
+    wEns(2) = nEl(1)/nEl(2)*wEns(2)
+    wEns(3) = nEl(1)/nEl(3)*wEns(3)
+
   end if
+
   write(*,*)'----------------------------------------------------------'
   write(*,*)' Ensemble weights '
   write(*,*)'----------------------------------------------------------'
@@ -123,8 +160,8 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
   
 ! Read parameters for weight-dependent functional
   read(1,*)
-  read(1,*) (aCC_w1(I),I=1,3)
-  read(1,*) (aCC_w2(I),I=1,3)
+  read(1,*) (aCC_w1(iParam),iParam=1,3)
+  read(1,*) (aCC_w2(iParam),iParam=1,3)
 ! Read choice of exchange coefficient
   read(1,*)
   read(1,*) Cx_choice
@@ -140,13 +177,6 @@ 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(*,*) 
- 
-! Read occupation numbers for orbitals nO and nO+1
-  read(1,*)
-  do J=1,nEns
-    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 cb5d4fd..9e4ffaf 100644
--- a/src/eDFT/restricted_auxiliary_energy.f90
+++ b/src/eDFT/restricted_auxiliary_energy.f90
@@ -12,7 +12,7 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
   integer,intent(in)            :: nEns
   integer,intent(in)            :: nO
   double precision,intent(in)   :: eps(nBas)
-  double precision,intent(in)   :: occnum(nspin,2,nEns)
+  double precision,intent(in)   :: occnum(nBas,nspin,nEns)
 
 
 ! Local variables
diff --git a/src/eDFT/restricted_density_matrix.f90 b/src/eDFT/restricted_density_matrix.f90
index 1911ab7..3c5eaba 100644
--- a/src/eDFT/restricted_density_matrix.f90
+++ b/src/eDFT/restricted_density_matrix.f90
@@ -12,7 +12,7 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P,occnum)
   integer,intent(in)            :: nEns
   integer,intent(in)            :: nO
   double precision,intent(in)   :: c(nBas,nBas)
-  double precision,intent(in)   :: occnum(nspin,2,nEns)
+  double precision,intent(in)   :: occnum(nBas,nspin,nEns)
 
 
 ! Local variables
diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90
index 7cda48b..eef0a41 100644
--- a/src/eDFT/restricted_individual_energy.f90
+++ b/src/eDFT/restricted_individual_energy.f90
@@ -38,7 +38,7 @@ 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)   :: occnum(nspin,2,nEns)
+  double precision,intent(in)   :: occnum(nBas,nspin,nEns)
   integer,intent(in)            :: Cx_choice
 
 
diff --git a/src/eDFT/unrestricted_auxiliary_energy.f90 b/src/eDFT/unrestricted_auxiliary_energy.f90
index 8d69e97..0ebec6d 100644
--- a/src/eDFT/unrestricted_auxiliary_energy.f90
+++ b/src/eDFT/unrestricted_auxiliary_energy.f90
@@ -12,85 +12,30 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
   integer,intent(in)            :: nEns
   integer,intent(in)            :: nO(nspin)
   double precision,intent(in)   :: eps(nBas,nspin)
-  double precision,intent(in)   :: occnum(nspin,2,nEns)
+  double precision,intent(in)   :: occnum(nBas,nspin,nEns)
 
 
 ! Local variables
 
-  integer                       :: iEns,ispin
+  integer                       :: iEns
+  integer                       :: ispin
+  integer                       :: p
 
 ! Output variables
 
   double precision,intent(out)  :: Eaux(nspin,nEns)
 
-!----------------------------------------------------------
-!-------------- GOK-UKS ----------------------------------
-!----------------------------------------------------------
-  iEns = 1
-  do ispin=1,nspin
-    Eaux(ispin,iEns) = sum(eps(1:nO(ispin),ispin))
+! Compute auxiliary energies for each state of the ensemble based on occupation numbers
+
+  Eaux(:,:) = 0d0
+  do iEns=1,nEns
+    do ispin=1,nspin
+      do p=1,nBas
+
+        Eaux(ispin,iEns) = Eaux(ispin,iEns) + occnum(p,ispin,iEns)*eps(p,ispin) 
+
+      end do
+    end do
   end do
 
-  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 
- 
-!  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 ab6c1bf..02f1d22 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,occnum)
+subroutine unrestricted_density_matrix(nBas,nEns,c,P,occnum)
 
 ! Calculate density matrices
 
@@ -10,89 +10,39 @@ subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P,occnum)
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nEns
-  integer,intent(in)            :: nO(nspin)
   double precision,intent(in)   :: c(nBas,nBas,nspin)
-  double precision,intent(in)   :: occnum(nspin,2,nEns)
+  double precision,intent(in)   :: occnum(nBas,nspin,nEns)
 
 
 ! Local variables
 
   integer                       :: ispin
   integer                       :: iEns
+  integer                       :: q
+  integer                       :: mu,nu
 
 ! Output variables
 
   double precision,intent(out)  :: P(nBas,nBas,nspin,nEns)
 
-!-------------------------------------------------------
-!------------------------ GOK-UKS ----------------------
-!-------------------------------------------------------
-  iEns = 1
-  do ispin=1,nspin
-    P(:,:,ispin,iEns) =  matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
+! Compute density matrix for each state of the ensemble based on occupation numbers
+
+  P(:,:,:,:) = 0d0
+  do iEns=1,nEns
+    do ispin=1,nspin
+      do mu=1,nBas
+        do nu=1,nBas
+          do q=1,nBas
+
+            P(mu,nu,ispin,iEns) = P(mu,nu,ispin,iEns) & 
+                                + occnum(q,ispin,iEns)*c(mu,q,ispin)*c(nu,q,ispin)
+
+          end do 
+        end do 
+      end do 
+    end do 
   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
-! (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)))
 
 end subroutine unrestricted_density_matrix
diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90
index f710040..228be81 100644
--- a/src/eDFT/unrestricted_individual_energy.f90
+++ b/src/eDFT/unrestricted_individual_energy.f90
@@ -41,7 +41,7 @@ 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)   :: occnum(nspin,2,nEns)
+  double precision,intent(in)   :: occnum(nBas,nspin,nEns)
   integer,intent(in)            :: Cx_choice