From a96719d0c2cbbba9fca029e3022badf97489b004 Mon Sep 17 00:00:00 2001
From: Clotilde Marut <clotildemarut@lcpqmc80.ups-tlse.fr>
Date: Tue, 29 Sep 2020 09:43:01 +0200
Subject: [PATCH] scaling of N-centered

---
 input/basis                                   | 44 +++++--------
 input/dft                                     | 26 ++++----
 input/molecule                                |  5 +-
 input/molecule.xyz                            |  5 +-
 input/weight                                  | 44 +++++--------
 scan_w.sh                                     | 50 ++++++++++-----
 ..._lda_exchange_derivative_discontinuity.f90 | 64 ++++++-------------
 .../UCC_lda_exchange_individual_energy.f90    | 56 ++++++----------
 .../US51_lda_exchange_individual_energy.f90   | 26 +++++++-
 src/eDFT/eDFT.f90                             |  9 +--
 src/eDFT/eDFT_UKS.f90                         |  6 +-
 src/eDFT/unrestricted_auxiliary_energy.f90    | 21 +++++-
 ...cted_exchange_derivative_discontinuity.f90 |  6 +-
 ...nrestricted_exchange_individual_energy.f90 |  5 +-
 src/eDFT/unrestricted_individual_energy.f90   | 49 +++++++++++---
 ..._lda_exchange_derivative_discontinuity.f90 |  9 ++-
 ...tricted_lda_exchange_individual_energy.f90 |  9 ++-
 17 files changed, 236 insertions(+), 198 deletions(-)

diff --git a/input/basis b/input/basis
index b2b2293..fb05e68 100644
--- a/input/basis
+++ b/input/basis
@@ -1,30 +1,18 @@
-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        
+1   3
+S   3
+  1     13.0100000              0.0196850        
+  2      1.9620000              0.1379770        
+  3      0.4446000              0.4781480        
 S   1
-  1      0.0280500              1.0000000        
-P   3
-  1      1.5340000              0.0227840        
-  2      0.2749000              0.1391070        
-  3      0.0736200              0.5003750        
+  1      0.1220000              1.0000000        
 P   1
-  1      0.0240300              1.0000000        
-D   1
-  1      0.1239000              1.0000000
-
+  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
diff --git a/input/dft b/input/dft
index 3d1826c..adfe387 100644
--- a/input/dft
+++ b/input/dft
@@ -18,22 +18,22 @@
   1
 # Number of states in ensemble (nEns)
   3
-# occupation numbers of orbitals
-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 
-
-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 
+# occupation numbers of orbitals nO and nO+1
+ 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 
+ 
+ 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 
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
-  0.00 1.00
+  0.00 1.00 
 # Ncentered ? 0 for NO 
- 1 
+ 0 
 # Parameters for CC weight-dependent exchange functional
-0.420431 0.069097 -0.295049
-0.135075 -0.00770826 -0.028057
+0.445525 0.0901503 -0.286898
+0.191734 -0.0364788 -0.017035
 # choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2
 2
 # GOK-DFT: maxSCF thresh   DIIS n_diis guess_type ortho_type
diff --git a/input/molecule b/input/molecule
index 058d6dd..81c624a 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,4 +1,5 @@
 # nAt nEla nElb nCore nRyd
-  1   2     1   0     0
+ 2   1     1   0     0
 # Znuc   x            y           z
-  Li     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 c9a5a65..6edc99d 100644
--- a/input/molecule.xyz
+++ b/input/molecule.xyz
@@ -1,3 +1,4 @@
-  1
+  2
 
- Li     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 b2b2293..fb05e68 100644
--- a/input/weight
+++ b/input/weight
@@ -1,30 +1,18 @@
-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        
+1   3
+S   3
+  1     13.0100000              0.0196850        
+  2      1.9620000              0.1379770        
+  3      0.4446000              0.4781480        
 S   1
-  1      0.0280500              1.0000000        
-P   3
-  1      1.5340000              0.0227840        
-  2      0.2749000              0.1391070        
-  3      0.0736200              0.5003750        
+  1      0.1220000              1.0000000        
 P   1
-  1      0.0240300              1.0000000        
-D   1
-  1      0.1239000              1.0000000
-
+  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
diff --git a/scan_w.sh b/scan_w.sh
index 799d780..e79e50a 100755
--- a/scan_w.sh
+++ b/scan_w.sh
@@ -12,17 +12,27 @@ w1=0.00
 XF=$3
 CF=$4
 
+# for H
+#aw1="1.49852 7.79815 25.1445"
+#aw2="0.424545 -0.0382349 -0.32472"
+
+
 # for He
-aw1="0.420431 0.069097 -0.295049"
-aw2="0.135075 -0.00770826 -0.028057"
+#aw1="0.429447 0.053506 -0.339391"
+#aw2="0.254939 -0.0893396 0.00765453"
 
 # for H2
-#aw1="0.445525 0.0901503 -0.286898"
-#aw2="0.191734 -0.0364788 -0.017035"
+aw1="0.445525 0.0901503 -0.286898"
+aw2="0.191734 -0.0364788 -0.017035"
 
 # for Li
-#aw1="0.0560976 -0.00796904 -0.0238162"
-#aw2="0.0360106 0.00979306 -0.0172286"
+#aw1="0.055105 -0.00943825 -0.0267771"
+#aw2="0.0359827 0.0096623 -0.0173542"
+
+# for Li+
+#aw1="0.503566, 0.137076, -0.348529"
+#aw2="0.0553828, 0.00830375, -0.0234602"
+
 
 # for B
 #aw1="0.052676 -0.00624118 -0.000368825"
@@ -36,6 +46,9 @@ aw2="0.135075 -0.00770826 -0.028057"
 #aw1="-0.00201219 -0.00371002 -0.00212719"
 #aw2="-0.00117715 0.00188738 -0.000414761"
 
+# for Be
+#aw1="0.0663282, -0.0117682, -0.0335909"
+#aw2="0.0479262, 0.00966351, -0.0208712"
 
 
 DATA=${MOL}_${BASIS}_${XF}_${CF}_${w2}.dat
@@ -65,6 +78,15 @@ do
   echo "  1" >> input/dft
   echo "# Number of states in ensemble (nEns)" >> input/dft
   echo "  3" >> input/dft
+  echo "# occupation numbers of orbitals nO and nO+1"  >> input/dft
+  echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
+  echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
+  echo " " >> input/dft
+  echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
+  echo " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
+  echo " " >> input/dft
+  echo " 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
+  echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
   echo "# Ensemble weights: wEns(1),...,wEns(nEns-1)" >> input/dft
   echo "  ${w1} ${w2} " >> input/dft
   echo "# Ncentered ? 0 for NO " >> input/dft
@@ -74,13 +96,6 @@ do
   echo ${aw2} >> input/dft
   echo "# choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2" >> input/dft
   echo "2"  >> input/dft
-  echo "# occupation numbers of orbitals nO and nO+1"  >> 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 " 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
@@ -93,7 +108,12 @@ do
   EA=`grep "Electronic Affinity"  ${OUTPUT} | grep " au" | tail -1 | cut -d":" -f 2 | sed 's/au//'`
   FG=`grep "Fundamental Gap"  ${OUTPUT} | grep " au" | tail -1 | cut -d":" -f 2 | sed 's/au//'`
   Ex=`grep "Exchange        energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/au//'` 
-  echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA $FG $Ex
-  echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA $FG $Ex >> ${DATA}
+  HOMOa=`grep "HOMO a    energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
+  LUMOa=`grep "LUMO a    energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
+  HOMOb=`grep "HOMO a    energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
+  LUMOb=`grep "LUMO b    energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
+
+  echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA $FG $Ex $HOMOa $LUMOa $HOMOb $LUMOb
+  echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA $FG $Ex $HOMOa $LUMOa $HOMOb $LUMOb >> ${DATA}
 done
 
diff --git a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90
index 1246d3b..6e164d8 100644
--- a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90
+++ b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90
@@ -1,4 +1,4 @@
-subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD,Cx_choice)
+subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,Cx_choice,doNcentered,ExDD)
 
 ! Compute the unrestricted version of the curvature-corrected exchange ensemble derivative
 
@@ -15,6 +15,8 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
   double precision,intent(in)   :: weight(nGrid)
   double precision,intent(in)   :: rhow(nGrid)
   integer,intent(in)            :: Cx_choice
+  integer,intent(in)            :: doNcentered
+
 
 ! Local variables
 
@@ -28,51 +30,39 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
   double precision              :: a2,b2,c2,w2
   double precision              :: dCxdw1,dCxdw2
 
+  double precision              :: nEli,nElw
+
+
 ! Output variables
 
   double precision,intent(out)  :: ExDD(nEns)
 
+! External variable
+
+  double precision,external     :: electron_number
+
+
 ! Memory allocation
 
   allocate(dExdw(nEns))
 
-! Single excitation parameters
 
-!  a1 = 0.0d0
-!  b1 = 0.0d0
-!  c1 = 0.0d0
-
-! Parameters for H2 at equilibrium
-
-! a2 = +0.5751782560799208d0
-! b2 = -0.021108186591137282d0
-! c2 = -0.36718902716347124d0
-
-! Parameters for stretch H2
-
-!  a2 = + 0.01922622507087411d0
-!  b2 = - 0.01799647558018601d0
-!  c2 = - 0.022945430666782573d0
-
-! Parameters for He
-
-! a2 = 1.9125735895875828d0
-! b2 = 2.715266992840757d0
-! c2 = 2.1634223380633086d0
-
-! Parameters for He N -> N-1
+! Parameters for N -> N-1
 
   a1 = aCC_w1(1)
   b1 = aCC_w1(2)
   c1 = aCC_w1(3)
 
 
-! Parameters for He N -> N+1
+! Parameters for N -> N+1
 
   a2 = aCC_w2(1)
   b2 = aCC_w2(2)
   c2 = aCC_w2(3)
  
+  
+  nElw = electron_number(nGrid,weight,rhow)
+
 
 ! Cx coefficient for unrestricted Slater LDA exchange
 
@@ -96,23 +86,6 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
            * (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))  
   end select
 
-! 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))
-
-!  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.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
@@ -143,4 +116,9 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
     end do
   end do
 
+  if (doNcentered .NE. 0) then
+    ExDD(2) = ((nElw-1)/nElw)*ExDD(2)
+    ExDD(3) = ((nElw+1)/nElw)*ExDD(3)
+  end if
+
 end subroutine UCC_lda_exchange_derivative_discontinuity
diff --git a/src/eDFT/UCC_lda_exchange_individual_energy.f90 b/src/eDFT/UCC_lda_exchange_individual_energy.f90
index 72867b3..019564d 100644
--- a/src/eDFT/UCC_lda_exchange_individual_energy.f90
+++ b/src/eDFT/UCC_lda_exchange_individual_energy.f90
@@ -1,4 +1,4 @@
-subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
+subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Cx_choice,doNcentered,Ex)
 
 ! Compute the unrestricted version of the curvature-corrected exchange functional
 
@@ -16,12 +16,14 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
   double precision,intent(in)   :: rhow(nGrid)
   double precision,intent(in)   :: rho(nGrid)
   integer,intent(in)            :: Cx_choice
+  integer,intent(in)            :: doNcentered
  
 ! Local variables
 
   integer                       :: iG
   double precision              :: r,rI,alpha
   double precision              :: e_p,dedr
+  double precision              :: nEli,nElw
 
   double precision              :: a1,b1,c1,w1
   double precision              :: a2,b2,c2,w2
@@ -31,37 +33,18 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
 
   double precision,intent(out)  :: Ex
 
-! Single excitation parameter
+! External variable
 
-!  a1 = 0.0d0
-!  b1 = 0.0d0
-!  c1 = 0.0d0
+  double precision,external     :: electron_number
 
-! Parameters for H2 at equilibrium
 
-! a2 = +0.5751782560799208d0
-! b2 = -0.021108186591137282d0
-! c2 = -0.36718902716347124d0
-
-! Parameters for stretch H2
-
-!  a2 = + 0.01922622507087411d0
-!  b2 = - 0.01799647558018601d0
-!  c2 = - 0.022945430666782573d0
-
-! Parameters for He
-
-! a2 = 1.9125735895875828d0
-! b2 = 2.715266992840757d0
-! c2 = 2.1634223380633086d0
-
-! Parameters for He N -> N-1
+! Parameters for N -> N-1
 
   a1 = aCC_w1(1)
   b1 = aCC_w1(2)
   c1 = aCC_w1(3)
 
-! Parameters for He N -> N+1
+! Parameters for N -> N+1
 
   a2 = aCC_w2(1)
   b2 = aCC_w2(2)
@@ -86,14 +69,10 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
   Cx = alpha*Fx2*Fx1
   end select
 
-! for two-weight ensembles
-!  Cx = alpha*Fx1*Fx2
+  nEli = electron_number(nGrid,weight,rho)
 
-! for left ensembles
-!  Cx = alpha*Fx1
-
-! for right ensembles
-!  Cx = alpha*Fx2  
+  nElw = electron_number(nGrid,weight,rhow)
+ 
 
 ! Compute LDA exchange matrix in the AO basis
 
@@ -107,13 +86,20 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
 
       e_p  =         Cx*r**(1d0/3d0)
       dedr = 1d0/3d0*Cx*r**(-2d0/3d0)
-
-      Ex = Ex - weight(iG)*dedr*r*r
+     
+      if (doNcentered == 0) then
+        Ex = Ex - weight(iG)*dedr*r*r
+      else
+        Ex = Ex - weight(iG)*dedr*r*r*(nEli/nElw)
+      end if
 
       if(rI > threshold) then
 
-      Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI)
-
+        if (doNcentered == 0) then
+          Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI)
+        else
+          Ex = Ex + weight(iG)*((nEli/nElw)*e_p*rI + dedr*r*rI)
+        end if
       endif
 
     endif
diff --git a/src/eDFT/US51_lda_exchange_individual_energy.f90 b/src/eDFT/US51_lda_exchange_individual_energy.f90
index 139be0f..389af78 100644
--- a/src/eDFT/US51_lda_exchange_individual_energy.f90
+++ b/src/eDFT/US51_lda_exchange_individual_energy.f90
@@ -1,4 +1,4 @@
-subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
+subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,Ex)
 
 ! Compute the restricted version of Slater's LDA exchange individual energy
 
@@ -11,17 +11,29 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
   double precision,intent(in)   :: weight(nGrid)
   double precision,intent(in)   :: rhow(nGrid)
   double precision,intent(in)   :: rho(nGrid)
+  integer,intent(in)            :: doNcentered
+
 
 ! Local variables
 
   integer                       :: iG
   double precision              :: r,rI,alpha
   double precision              :: e,dedr
+  double precision              :: nEli,nElw
 
 ! Output variables
 
   double precision,intent(out)  :: Ex
 
+! External variable
+
+  double precision,external     :: electron_number
+
+  nEli = electron_number(nGrid,weight,rho)
+
+  nElw = electron_number(nGrid,weight,rhow)
+
+
 ! Compute LDA exchange matrix in the AO basis
 
   alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
@@ -37,11 +49,19 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
       e    =         alpha*r**(1d0/3d0)
       dedr = 1d0/3d0*alpha*r**(-2d0/3d0)
 
-      Ex = Ex - weight(iG)*dedr*r*r
+      if (doNcentered == 0) then
+        Ex = Ex - weight(iG)*dedr*r*r      
+      else
+        Ex = Ex - weight(iG)*dedr*r*r*(nEli/nElw)
+      end if
 
       if(rI > threshold) then
 
-      Ex = Ex + weight(iG)*(e*rI + dedr*r*rI)
+        if (doNcentered == 0) then
+          Ex = Ex + weight(iG)*(e*rI + dedr*r*rI)
+        else
+          Ex = Ex + weight(iG)*((nEli/nElw)*e*rI + dedr*r*rI)
+        end if
 
       endif
 
diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90
index f4e2f29..092bdab 100644
--- a/src/eDFT/eDFT.f90
+++ b/src/eDFT/eDFT.f90
@@ -187,7 +187,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(:,:),occnum,Cx_choice)
+                 S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered)
     call cpu_time(end_KS)
 
     t_KS = end_KS - start_KS
@@ -205,7 +205,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(:,:),occnum,Cx_choice)
+                 S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered)
     call cpu_time(end_KS)
 
     t_KS = end_KS - start_KS
@@ -222,7 +222,8 @@ 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,occnum,Cx_choice)
+                   nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum, &
+                   Cx_choice,doNcentered)
     call cpu_time(end_KS)
 
     t_KS = end_KS - start_KS
@@ -239,7 +240,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,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice)
+                  nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice,doNcentered)
     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 53b5b99..d5e3d2a 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,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice)
+                   nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice,doNcentered)
 
 ! Perform unrestricted Kohn-Sham calculation for ensembles
 
@@ -31,6 +31,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
   double precision,intent(in)   :: ENuc
   double precision,intent(in)   :: occnum(nBas,nspin,nEns)
   integer,intent(in)            :: Cx_choice
+  integer,intent(in)            :: doNcentered
 
 ! Local variables
 
@@ -380,6 +381,7 @@ 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,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum,Cx_choice)
+                                      AO,dAO,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum, &
+                                      Cx_choice,doNcentered)
 
 end subroutine eDFT_UKS
diff --git a/src/eDFT/unrestricted_auxiliary_energy.f90 b/src/eDFT/unrestricted_auxiliary_energy.f90
index d493d2b..34a6e8e 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,eps,Eaux,occnum)
+subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,doNcentered,Eaux)
 
 ! Compute the auxiliary KS energies 
 
@@ -12,18 +12,28 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,Eaux,occnum)
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: eps(nBas,nspin)
   double precision,intent(in)   :: occnum(nBas,nspin,nEns)
-
+  integer,intent(in)            :: doNcentered
 
 ! Local variables
 
-  integer                       :: iEns
+  integer                       :: iEns,iBas
   integer                       :: ispin
   integer                       :: p
+  double precision,allocatable  :: nEl(:) 
+
 
 ! Output variables
 
   double precision,intent(out)  :: Eaux(nspin,nEns)
 
+  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
+
 ! Compute auxiliary energies for each state of the ensemble based on occupation numbers
 
   Eaux(:,:) = 0d0
@@ -35,6 +45,11 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,Eaux,occnum)
 
       end do
     end do
+
+!    if (doNcentered .NE. 0) then
+!      Eaux(:,iEns) = Eaux(:,iEns)*(nEl(iEns)/nEl(1))   
+!    end if 
+
   end do
 
 end subroutine unrestricted_auxiliary_energy
diff --git a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90
index 30143bc..39547c7 100644
--- a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90
+++ b/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90
@@ -1,4 +1,5 @@
-subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,ExDD,Cx_choice)
+subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,&
+                                                          Cx_choice,doNcentered,ExDD)
 
 ! Compute the exchange part of the derivative discontinuity
 
@@ -18,6 +19,7 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC
   double precision,intent(in)   :: rhow(nGrid)
   double precision,intent(in)   :: drhow(ncart,nGrid)
   integer,intent(in)            :: Cx_choice
+  integer,intent(in)            :: doNcentered
 
 ! Local variables
 
@@ -39,7 +41,7 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC
     case(1) 
 
       call unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),&
-                                                              rhow(:),ExDD(:),Cx_choice)
+                                                              rhow(:),Cx_choice,doNcentered,ExDD(:))
 
 !   GGA functionals
 
diff --git a/src/eDFT/unrestricted_exchange_individual_energy.f90 b/src/eDFT/unrestricted_exchange_individual_energy.f90
index b0d5a66..dba469f 100644
--- a/src/eDFT/unrestricted_exchange_individual_energy.f90
+++ b/src/eDFT/unrestricted_exchange_individual_energy.f90
@@ -1,5 +1,5 @@
 subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, & 
-                                                   ERI,Pw,P,rhow,drhow,rho,drho,Ex,Cx_choice)
+                                                   ERI,Pw,P,rhow,drhow,rho,drho,Cx_choice,doNcentered,Ex)
 
 ! Compute the exchange individual energy
 
@@ -26,6 +26,7 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE
   double precision,intent(in)   :: rho(nGrid)
   double precision,intent(in)   :: drho(ncart,nGrid)
   integer,intent(in)            :: Cx_choice
+  integer,intent(in)            :: doNcentered
 
 ! Local variables
 
@@ -50,7 +51,7 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE
     case(1) 
 
       call unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,&
-                                                       rhow,rho,ExLDA,Cx_choice)
+                                                       rhow,rho,Cx_choice,doNcentered,ExLDA)
 
       Ex = ExLDA
 
diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90
index 2222126..10efd57 100644
--- a/src/eDFT/unrestricted_individual_energy.f90
+++ b/src/eDFT/unrestricted_individual_energy.f90
@@ -1,5 +1,6 @@
 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, &
-                                          T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum,Cx_choice)
+                                          T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum,&
+                                          Cx_choice,doNcentered)
 
 ! Compute unrestricted individual energies as well as excitation energies
 
@@ -42,6 +43,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
   double precision              :: Ew
   double precision,intent(in)   :: occnum(nBas,nspin,nEns)
   integer,intent(in)            :: Cx_choice
+  integer,intent(in)            :: doNcentered
 
 
 ! Local variables
@@ -64,7 +66,9 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
 
   double precision,external     :: trace_matrix
 
-  integer                       :: ispin,iEns
+  integer                       :: ispin,iEns,iBas
+  double precision,allocatable  :: nEl(:)
+
   
   double precision,external     :: electron_number
 
@@ -73,15 +77,31 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
   double precision,intent(out)  :: E(nEns)
   double precision,intent(out)  :: Om(nEns)
 
+ 
+  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*,'test1'
+
 !------------------------------------------------------------------------
 ! Kinetic energy
 !------------------------------------------------------------------------
 
   do ispin=1,nspin
     do iEns=1,nEns
-      ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:)))
+      if (doNcentered == 0) then
+        ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) 
+      else 
+        ET(ispin,iEns) = (nEl(iEns)/nEl(1))*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:)))
+      end if
     end do
   end do
+  print*,'test2'
 
 !------------------------------------------------------------------------
 ! Potential energy
@@ -89,10 +109,15 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
 
   do iEns=1,nEns
     do ispin=1,nspin
-      EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
-    end do
+      if (doNcentered == 0) then
+        EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
+      else 
+        EV(ispin,iEns) = (nEl(iEns)/nEl(1))*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
+      end if  
+  end do
   end do
 
+  print*,'test3'
 !------------------------------------------------------------------------
 ! Individual Hartree energy
 !------------------------------------------------------------------------
@@ -113,9 +138,11 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
 
     EJ(3,iEns) =       trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) &
                - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
-
+    if (doNcentered .NE. 0) then
+      EJ(:,iEns) = (nEl(iEns)/nEl(1))*EJ(:,iEns)
+    end if
   end do
-     
+       print*,'test4'
 !------------------------------------------------------------------------
 ! Checking Hartree contributions for each individual states
 !------------------------------------------------------------------------
@@ -139,10 +166,12 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
     do ispin=1,nspin
       call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,ERI, &
                                                    Pw(:,:,ispin),P(:,:,ispin,iEns),rhow(:,ispin),drhow(:,:,ispin),          &  
-                                                   rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns),Cx_choice)
+                                                   rho(:,ispin,iEns),drho(:,:,ispin,iEns),Cx_choice,doNcentered,Ex(ispin,iEns))
     end do
   end do
 
+  print*,'test5'
+
 !------------------------------------------------------------------------
 ! Checking exchange contributions for each individual states
 !------------------------------------------------------------------------
@@ -184,7 +213,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
 ! Compute auxiliary energies
 !------------------------------------------------------------------------
  
-  call unrestricted_auxiliary_energy(nBas,nEns,eps,Eaux,occnum)
+  call unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,doNcentered,Eaux)
 
 !------------------------------------------------------------------------
 ! Compute derivative discontinuities
@@ -193,7 +222,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
   do ispin=1,nspin 
 
     call unrestricted_exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, &
-                                                        rhow(:,ispin),drhow(:,:,ispin),ExDD(ispin,:),Cx_choice)
+                                                        rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered,ExDD(ispin,:))
   end do
 
   call unrestricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD)
diff --git a/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90
index 4750bd1..0e80f3a 100644
--- a/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90
+++ b/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90
@@ -1,4 +1,5 @@
-subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD,Cx_choice)
+subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,&
+                                                              Cx_choice,doNcentered,ExDD)
 
 ! Compute the exchange LDA part of the derivative discontinuity
 
@@ -17,6 +18,7 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_
   double precision,intent(in)   :: weight(nGrid)
   double precision,intent(in)   :: rhow(nGrid)
   integer,intent(in)            :: Cx_choice
+  integer,intent(in)            :: doNcentered
 
 ! Local variables
 
@@ -25,7 +27,7 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_
 
   double precision,intent(out)  :: ExDD(nEns)
 
-! Select correlation functional
+! Select exchange functional
 
   select case (DFA)
 
@@ -35,7 +37,8 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_
 
     case ('CC')
 
-      call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:),Cx_choice)
+      call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),&
+                                                     Cx_choice,doNcentered,ExDD(:))
 
     case default
 
diff --git a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90
index 9d64dfe..9ac05b9 100644
--- a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90
+++ b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90
@@ -1,4 +1,5 @@
-subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
+subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,&
+                                                       Cx_choice,doNcentered,Ex)
 
 ! Compute LDA exchange energy for individual states
 
@@ -18,6 +19,7 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn
   double precision,intent(in)   :: rhow(nGrid)
   double precision,intent(in)   :: rho(nGrid)
   integer,intent(in)            :: Cx_choice
+  integer,intent(in)            ::doNcentered
 
 ! Output variables
 
@@ -29,11 +31,12 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn
 
     case ('S51')
 
-      call US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
+      call US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,Ex)
 
     case ('CC')
 
-      call UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
+      call UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,&
+                                              Cx_choice,doNcentered,Ex)
 
     case default