diff --git a/include/parameters.h b/include/parameters.h
index 89c2c93..2644d6a 100644
--- a/include/parameters.h
+++ b/include/parameters.h
@@ -2,6 +2,7 @@
   integer,parameter              :: nspin = 2
   integer,parameter              :: nsp = 3
   integer,parameter              :: maxEns = 4
+  integer,parameter              :: maxCC = 5
   integer,parameter              :: maxShell = 512
   integer,parameter              :: maxL = 7
   integer,parameter              :: n1eInt = 3
diff --git a/input/dft b/input/dft
index 2be05f9..85cd470 100644
--- a/input/dft
+++ b/input/dft
@@ -4,38 +4,38 @@
 # Hartree      = 0: H
 # LDA          = 1: S51,CC-S51
 # GGA          = 2: B88,G96,PBE
-# MGGA         = 3:
-# Hybrid       = 4 HF,B3LYP,PBE
-1 S51  
+# MGGA         = 3: 
+# Hybrid       = 4: HF,B3LYP,PBE
+  1 S51
 # correlation rung: 
 # Hartree      = 0: H
 # LDA          = 1: PW92,VWN3,VWN5,eVWN5
 # GGA          = 2: LYP,PBE
 # MGGA         = 3:
-# Hybrid       = 4: HF,B3LYP,PBE 
-1 VWN5 
+# Hybrid       = 4: HF,B3LYP,PBE
+ 1 VWN5
 # quadrature grid SG-n
-  1
+  0
 # Number of states in ensemble (nEns)
-2
-# occupation numbers
-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 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 0 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 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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+  2
+# occupation numbers 
+ 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 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 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 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 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 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 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 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 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 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 0 0 0 0 0 0 0 0 0 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 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 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 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 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 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 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 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 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 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 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 0 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 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 0 0 0 0 0 0 0 0
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
-  0.95 0.0 0.0
-# Ncentered ? 
-F
+ 0.95 0.00 0.00
+# N-centered?
+ F   
 # Parameters for CC weight-dependent exchange functional
--0.766201 -0.155585 0.00130104 0.0
-0.00 0.00 0.00 0.00 0.0
+3
+0.0 0.0 0.0
 # choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2
-1
+2
diff --git a/input/options b/input/options
index f54eee2..2592c87 100644
--- a/input/options
+++ b/input/options
@@ -3,16 +3,16 @@
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis
-      64     0.0000000001  T    5
+      64     0.0000001  T    5
 # spin: TDA singlet triplet spin_conserved spin_flip 
          F   T       T       T              T 
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm 
        256  0.00001 T    5      T    0.00367493  3      
 # GW/GT: maxSCF thresh  DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 
-         256   0.00001   T    5     T   0.00367493 F      F     F     F   F  
+         256   0.00001   T    5     T   0.00 F      F     F     F   F  
 # ACFDT: AC Kx  XBS
-         T  T   F
+         F  T   F
 # BSE: BSE dBSE dTDA evDyn
-        T    F    T    F
+        T    T    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T
diff --git a/mol/h2.xyz b/mol/h2.xyz
index a4e936a..83cc4a4 100644
--- a/mol/h2.xyz
+++ b/mol/h2.xyz
@@ -1,4 +1,4 @@
 2
  
 H 0. 0. 0.
-H 0. 0. 2.000000
+H 0. 0. 3.2
diff --git a/src/MBPT/G0T0.f90 b/src/MBPT/G0T0.f90
index a838c40..3054daa 100644
--- a/src/MBPT/G0T0.f90
+++ b/src/MBPT/G0T0.f90
@@ -198,6 +198,9 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
 
   if(BSE) then
 
+!   eG0T0(1) = -0.5507952119d0
+!   eG0T0(2) = +1.540259769d0
+
     call Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,   & 
                                 Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t, &
                                 ERI_MO,dipole_int,eHF,eG0T0,EcBSE)
diff --git a/src/MBPT/dynamic_Tmatrix_A.f90 b/src/MBPT/dynamic_Tmatrix_A.f90
index 55d4c44..8472e0f 100644
--- a/src/MBPT/dynamic_Tmatrix_A.f90
+++ b/src/MBPT/dynamic_Tmatrix_A.f90
@@ -62,10 +62,10 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
           end do
 
             do kl=1,nOO
-            chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
+            chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
           end do
 
-          A_dyn(ia,jb) = A_dyn(ia,jb) - 1d0*lambda*chi
+          A_dyn(ia,jb) = A_dyn(ia,jb) + 1d0*lambda*chi
 
           chi = 0d0
 
@@ -75,11 +75,11 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
           end do
 
           do kl=1,nOO
-            eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
-            chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
+            eps = + OmBSE - Omega2(kl) - (eGT(a) + eGT(b))
+            chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
           end do
 
-          A_dyn(ia,jb) = A_dyn(ia,jb) + 1d0*lambda*chi
+          A_dyn(ia,jb) = A_dyn(ia,jb) - 1d0*lambda*chi
 
           chi = 0d0
 
@@ -89,11 +89,11 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
           end do
 
           do kl=1,nOO
-            eps = + OmBSE - Omega2(kl) - (eGT(a) + eGT(b))
+            eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
             chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
           end do
 
-          ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - 1d0*lambda*chi
+          ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + 1d0*lambda*chi
 
         end do
       end do
diff --git a/src/MBPT/static_Tmatrix_TA.f90 b/src/MBPT/static_Tmatrix_TA.f90
index c3f2583..1c78e0d 100644
--- a/src/MBPT/static_Tmatrix_TA.f90
+++ b/src/MBPT/static_Tmatrix_TA.f90
@@ -49,11 +49,11 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
           enddo
 
           do kl=1,nOO
-!           chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
+!           chi = chi + lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
             chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
           enddo
 
-          TA(ia,jb) = TA(ia,jb) + 1d0*lambda*chi
+          TA(ia,jb) = TA(ia,jb) - 1d0*lambda*chi
 
         enddo
       enddo
diff --git a/src/MBPT/static_Tmatrix_TB.f90 b/src/MBPT/static_Tmatrix_TB.f90
index d502c89..d1f993c 100644
--- a/src/MBPT/static_Tmatrix_TB.f90
+++ b/src/MBPT/static_Tmatrix_TB.f90
@@ -53,7 +53,7 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
             chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
           enddo
 
-          TB(ia,jb) = TB(ia,jb) + 1d0*lambda*chi
+          TB(ia,jb) = TB(ia,jb) - 1d0*lambda*chi
 
         enddo
       enddo
diff --git a/src/RPA/ACFDT_Tmatrix.f90 b/src/RPA/ACFDT_Tmatrix.f90
index 3245e34..b667c15 100644
--- a/src/RPA/ACFDT_Tmatrix.f90
+++ b/src/RPA/ACFDT_Tmatrix.f90
@@ -64,11 +64,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
 
 ! Useful quantities
 
-! nOOs = nO*nO
-! nVVs = nV*nV
-
-  nOOs = nO*(nO+1)/2
-  nVVs = nV*(nV+1)/2
+  nOOs = nO*nO
+  nVVs = nV*nV
 
   nOOt = nO*(nO-1)/2
   nVVt = nV*(nV-1)/2
@@ -121,31 +118,31 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
       TA(:,:) = 0d0
       TB(:,:) = 0d0
 
-!     if(doXBS) then
+      if(doXBS) then
 
-!       isp_T  = 1
-!       iblock = 3
+        isp_T  = 1
+        iblock = 3
 
-!       call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI,  &
-!                               Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T))
+        call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI,  &
+                                Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T))
 
-!       call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
+        call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
-!                    call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
-!       if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
+                     call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
+        if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
 
-!       isp_T  = 2
-!       iblock = 4
+        isp_T  = 2
+        iblock = 4
 
-!       call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI,  &
-!                               Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T))
+        call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI,  &
+                                Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T))
 
-!       call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
+        call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
-!                    call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
-!       if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
+                     call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
+        if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
 
-!     end if
+      end if
 
       call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, &
                                    EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
@@ -191,31 +188,31 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
       TA(:,:) = 0d0
       TB(:,:) = 0d0
 
-!     if(doXBS) then
+      if(doXBS) then
 
-!       isp_T  = 1
-!       iblock = 3
+        isp_T  = 1
+        iblock = 3
 
-!       call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI,  &
-!                               Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T))
+        call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI,  &
+                                Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T))
 
-!       call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
+        call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
-!                    call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
-!       if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
+                     call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
+        if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
 
-!       isp_T  = 2
-!       iblock = 4
+        isp_T  = 2
+        iblock = 4
 
-!       call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI,  &
-!                               Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T))
+        call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI,  &
+                                Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T))
 
-!       call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
+        call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
-!                    call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
-!       if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
+                     call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
+        if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
 
-!     end if  
+      end if  
 
       call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, &
                                    EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
diff --git a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90
index 069c4b1..a8b0d1a 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,Cx_choice,doNcentered,kappa,ExDD)
+subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,doNcentered,kappa,ExDD)
 
 ! Compute the unrestricted version of the curvature-corrected exchange ensemble derivative
 
@@ -9,8 +9,8 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
 
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: wEns(nEns)
-  double precision,intent(in)   :: aCC_w1(3)
-  double precision,intent(in)   :: aCC_w2(3)
+  integer,intent(in)            :: nCC
+  double precision,intent(in)   :: aCC(nCC,nEns-1)
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   double precision,intent(in)   :: rhow(nGrid)
@@ -44,22 +44,22 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
   allocate(dExdw(nEns))
 
 
-! Parameters for N -> N-1
+! Parameters for first state
 
-  a1 = aCC_w1(1)
-  b1 = aCC_w1(2)
-  c1 = aCC_w1(3)
+  a1 = aCC(1,1)
+  b1 = aCC(2,1)
+  c1 = aCC(3,1)
 
+! Parameters for second state
 
-! Parameters for N -> N+1
-
-  a2 = aCC_w2(1)
-  b2 = aCC_w2(2)
-  c2 = aCC_w2(3)
+  a2 = aCC(1,2)
+  b2 = aCC(2,2)
+  c2 = aCC(3,2)
  
   w1 = wEns(2)
   w2 = wEns(3)
 
+! Defining enhancements factor for weight-dependent functionals
 
   if (doNcentered) then
   
@@ -88,7 +88,6 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
 
   else
 
-
     select case (Cx_choice)
 
       case(1)
diff --git a/src/eDFT/UCC_lda_exchange_energy.f90 b/src/eDFT/UCC_lda_exchange_energy.f90
index 51238ff..5bc917d 100644
--- a/src/eDFT/UCC_lda_exchange_energy.f90
+++ b/src/eDFT/UCC_lda_exchange_energy.f90
@@ -1,4 +1,4 @@
-subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Cx_choice,doNcentered,Ex)
+subroutine UCC_lda_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex)
 
 ! Compute the unrestricted version of the curvature-corrected exchange functional
 
@@ -9,8 +9,8 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Cx_c
 
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: wEns(nEns)
-  double precision,intent(in)   :: aCC_w1(3)
-  double precision,intent(in)   :: aCC_w2(3)
+  integer,intent(in)            :: nCC
+  double precision,intent(in)   :: aCC(nCC,nEns-1)
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   double precision,intent(in)   :: rho(nGrid)
@@ -30,44 +30,19 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Cx_c
 
   double precision              :: Ex
 
-! Single excitation parameter
+! Parameters for first state
 
-!  a1 = 0.0d0
-!  b1 = 0.0d0
-!  c1 = 0.0d0
+  a1 = aCC(1,1)
+  b1 = aCC(2,1)
+  c1 = aCC(3,1)
 
-! Parameters for H2 at equilibrium
+! Parameters for second state
 
-! a2 = +0.5751782560799208d0
-! b2 = -0.021108186591137282d0
-! c2 = -0.36718902716347124d0
+  a2 = aCC(1,2)
+  b2 = aCC(2,2)
+  c2 = aCC(3,2)
 
-! 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
-
-  a1 = aCC_w1(1)
-  b1 = aCC_w1(2)
-  c1 = aCC_w1(3)
-
-! Parameters for He N -> N+1
-
-  a2 = aCC_w2(1)
-  b2 = aCC_w2(2)
-  c2 = aCC_w2(3)
-
-! Fx1 for states N and N-1
-! Fx2 for states N and N+1
+! Defining enhancements factor for weight-dependent functionals
 
   if(doNcentered) then 
     
@@ -79,7 +54,6 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Cx_c
 
   else  
 
-
     w1 = wEns(2)
     Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)
 
@@ -112,9 +86,7 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Cx_c
 
     r = max(0d0,rho(iG))
 
-    if(r > threshold) then
-      Ex = Ex + weight(iG)*Cx*r**(4d0/3d0) 
-    endif
+    if(r > threshold) Ex = Ex + weight(iG)*Cx*r**(4d0/3d0) 
 
   enddo
 
diff --git a/src/eDFT/UCC_lda_exchange_individual_energy.f90 b/src/eDFT/UCC_lda_exchange_individual_energy.f90
index 3671707..cd0666e 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,Cx_choice,doNcentered,kappa,Ex)
+subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho,Cx_choice,doNcentered,kappa,Ex)
 
 ! Compute the unrestricted version of the curvature-corrected exchange functional
 
@@ -9,8 +9,8 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
 
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: wEns(nEns)
-  double precision,intent(in)   :: aCC_w1(3)
-  double precision,intent(in)   :: aCC_w2(3)
+  integer,intent(in)            :: nCC
+  double precision,intent(in)   :: aCC(nCC,nEns-1)
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   double precision,intent(in)   :: rhow(nGrid)
@@ -39,17 +39,19 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
   double precision,external     :: electron_number
 
 
-! Parameters for N -> N-1
+! Parameters for first state
 
-  a1 = aCC_w1(1)
-  b1 = aCC_w1(2)
-  c1 = aCC_w1(3)
+  a1 = aCC(1,1)
+  b1 = aCC(2,1)
+  c1 = aCC(3,1)
 
-! Parameters for N -> N+1
+! Parameters for second state
 
-  a2 = aCC_w2(1)
-  b2 = aCC_w2(2)
-  c2 = aCC_w2(3)
+  a2 = aCC(1,2)
+  b2 = aCC(2,2)
+  c2 = aCC(3,2)
+
+! Defining enhancements factor for weight-dependent functionals
 
   if(doNcentered) then
 
@@ -61,7 +63,6 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
 
   else
 
-
     w1 = wEns(2)
     Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)
 
diff --git a/src/eDFT/UCC_lda_exchange_potential.f90 b/src/eDFT/UCC_lda_exchange_potential.f90
index b8b12d5..114a204 100644
--- a/src/eDFT/UCC_lda_exchange_potential.f90
+++ b/src/eDFT/UCC_lda_exchange_potential.f90
@@ -1,4 +1,4 @@
-subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice,doNcentered)
+subroutine UCC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Fx,Cx_choice,doNcentered)
 
 ! Compute the unrestricted version of the curvature-corrected exchange potential
 
@@ -9,8 +9,8 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
 
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: wEns(nEns)
-  double precision,intent(in)   :: aCC_w1(3)
-  double precision,intent(in)   :: aCC_w2(3)
+  integer,intent(in)            :: nCC
+  double precision,intent(in)   :: aCC(nCC,nEns-1)
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   integer,intent(in)            :: nBas
@@ -32,44 +32,19 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
 
   double precision,intent(out)  :: Fx(nBas,nBas)
 
-! Single excitation parameter
+! Parameters for first state
 
-!  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
-
-  a1 = aCC_w1(1)
-  b1 = aCC_w1(2)
-  c1 = aCC_w1(3)
+  a1 = aCC(1,1)
+  b1 = aCC(2,1)
+  c1 = aCC(3,1)
  
-! Parameters for He N -> N+1
+! Parameters for second state
 
-  a2 = aCC_w2(1)
-  b2 = aCC_w2(2)
-  c2 = aCC_w2(3)
+  a2 = aCC(1,2)
+  b2 = aCC(2,2)
+  c2 = aCC(3,2)
 
-! Fx1 for states N and N-1
-! Fx2 for states N and N+1
+! Defining enhancements factor for weight-dependent functionals
 
   if(doNcentered) then 
 
diff --git a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90
index 19e4c9a..f26936e 100644
--- a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90
+++ b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90
@@ -23,7 +23,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
   double precision              :: a_f,x0_f,xx0_f,b_f,c_f,x_f,q_f
   double precision              :: a_a,x0_a,xx0_a,b_a,c_a,x_a,q_a
   double precision              :: dfzdz,dxdrs,dxdx_p,dxdx_f,dxdx_a,decdx_p,decdx_f,decdx_a
-  double precision              :: dzdr ,dfzdr ,drsdr ,decdr_p ,decdr_f ,decdr_a, decdr
+  double precision              :: dzdra,dzdrb,dfzdra,dfzdrb,drsdr,decdr_p,decdr_f,decdr_a,decdra,decdrb,decdr
   double precision              :: ec_z,ec_p,ec_f,ec_a
   double precision              :: fz,d2fz
 
@@ -139,11 +139,15 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
       ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
                  - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
    
-      ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+      ec_z = ec_p + ec_a*fz/d2fz*(1d0 - z**4) + (ec_f - ec_p)*fz*z**4
 
-      dzdr  = (1d0 - z)/r
       dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
-      dfzdr = dzdr*dfzdz
+
+      dzdra = + (1d0 - z)/r
+      dfzdra = dzdra*dfzdz
+
+      dzdrb = - (1d0 + z)/r
+      dfzdrb = dzdrb*dfzdz
 
       drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
       dxdrs = 0.5d0/sqrt(rs)
@@ -153,24 +157,29 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
       dxdx_a = 2d0*x + b_a
 
       decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
-                      - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+                      - b_p*x0_p/xx0_p*( 2d0/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
 
       decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
-                      - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+                      - b_f*x0_f/xx0_f*( 2d0/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
 
       decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
-                      - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+                      - b_a*x0_a/xx0_a*( 2d0/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
 
       decdr_p = drsdr*dxdrs*decdx_p
       decdr_f = drsdr*dxdrs*decdx_f
       decdr_a = drsdr*dxdrs*decdx_a
 
-      decdr = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdr/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdr*z**3 &
-            + (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdr*z**4 + 4d0*(ec_f - ec_p)*fz*dzdr*z**3
+      decdra = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 &
+             + (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
 
-      decdr  = ec_z + decdr*r
+      decdrb = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 &
+             + (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3
+      
+      decdr = 0d0
+      if(ra > threshold) decdr = decdr + decdra
+      if(rb > threshold) decdr = decdr + decdrb
  
-      Ecrr(2)  = Ecrr(2)  - weight(iG)*decdr*r*r
+      Ecrr(2)  = Ecrr(2) - weight(iG)*decdr*r*r
 
       if(rI > threshold) then
 
diff --git a/src/eDFT/UVWN5_lda_correlation_potential.f90 b/src/eDFT/UVWN5_lda_correlation_potential.f90
index c40ea59..0ffb8c0 100644
--- a/src/eDFT/UVWN5_lda_correlation_potential.f90
+++ b/src/eDFT/UVWN5_lda_correlation_potential.f90
@@ -34,7 +34,7 @@ subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
 
 ! Parameters of the functional
 
-  a_p  = +0.0621814D0/2D0
+  a_p  = +0.0621814D0/2d0
   x0_p = -0.10498d0
   b_p  = +3.72744d0
   c_p  = +12.9352d0
@@ -97,7 +97,7 @@ subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
        
           ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
 
-          dzdra = (1d0 - z)/r
+          dzdra = + (1d0 - z)/r
           dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
           dfzdra = dzdra*dfzdz
 
@@ -109,13 +109,13 @@ subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
           dxdx_a = 2d0*x + b_a
 
           decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
-                          - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+                          - b_p*x0_p/xx0_p*( 2d0/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
 
           decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
-                          - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+                          - b_f*x0_f/xx0_f*( 2d0/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
 
           decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
-                          - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+                          - b_a*x0_a/xx0_a*( 2d0/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
 
           decdra_p = drsdra*dxdrs*decdx_p
           decdra_f = drsdra*dxdrs*decdx_f
@@ -167,13 +167,13 @@ subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
           dxdx_a = 2d0*x + b_a
 
           decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
-                          - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+                          - b_p*x0_p/xx0_p*( 2d0/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
 
           decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
-                          - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+                          - b_f*x0_f/xx0_f*( 2d0/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
 
           decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
-                          - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+                          - b_a*x0_a/xx0_a*( 2d0/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
 
           decdrb_p = drsdrb*dxdrs*decdx_p
           decdrb_f = drsdrb*dxdrs*decdx_f
@@ -181,6 +181,7 @@ subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
 
           decdrb = decdrb_p + decdrb_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 &
                  + (decdrb_f - decdrb_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3
+
           Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdrb*r)
        
         end if
diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90
index 4133c1d..b58fa7f 100644
--- a/src/eDFT/eDFT.f90
+++ b/src/eDFT/eDFT.f90
@@ -62,8 +62,9 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
   integer                       :: nGrid
   double precision,allocatable  :: root(:,:)
   double precision,allocatable  :: weight(:)
-  double precision              :: aCC_w1(3)
-  double precision              :: aCC_w2(3)
+
+  integer                       :: nCC
+  double precision,allocatable  :: aCC(:,:)
 
   double precision,allocatable  :: AO(:,:)
   double precision,allocatable  :: dAO(:,:,:)
@@ -99,21 +100,14 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
   write(*,*) '******************************************'
   write(*,*)
 
-! Libxc version
-
-! call xc_f90_version(vmajor, vminor, vmicro)
-! write(*,'("Libxc version: ",I1,".",I1,".",I1)') vmajor, vminor, vmicro
-
-! call xcinfo()
-
 !------------------------------------------------------------------------
 ! DFT options
 !------------------------------------------------------------------------
 
 ! Allocate ensemble weights and MO coefficients
 
-  allocate(wEns(maxEns),occnum(nBas,nspin,maxEns))
-  call read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, & 
+  allocate(wEns(maxEns),aCC(maxCC,nEns-1),occnum(nBas,nspin,maxEns))
+  call read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,nCC,aCC, & 
                         doNcentered,occnum,Cx_choice)
 
 !------------------------------------------------------------------------
@@ -245,7 +239,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
     end do
 
     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,mix, & 
+    call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, & 
                   nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
     call cpu_time(end_KS)
 
@@ -262,7 +256,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
   if(method == 'eDFT-UKS') then
 
     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,mix, & 
+    call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, & 
                   nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
     call cpu_time(end_KS)
 
diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90
index 5326336..01f9217 100644
--- a/src/eDFT/eDFT_UKS.f90
+++ b/src/eDFT/eDFT_UKS.f90
@@ -1,4 +1,4 @@
-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,mix, &
+subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
                    nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eps,c,Pw,Vxc)
 
 ! Perform unrestricted Kohn-Sham calculation for ensembles
@@ -12,8 +12,8 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
   integer,intent(in)            :: x_DFA,c_DFA
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: wEns(nEns)
-  double precision,intent(in)   :: aCC_w1(3)
-  double precision,intent(in)   :: aCC_w2(3)
+  integer,intent(in)            :: nCC
+  double precision,intent(in)   :: aCC(nCC,nEns-1)
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   integer,intent(in)            :: maxSCF,max_diis,guess_type
@@ -259,7 +259,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
 !   Compute exchange potential
 
     do ispin=1,nspin
-      call unrestricted_exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,    &
+      call unrestricted_exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),nCC,aCC,nGrid,weight(:),nBas,    &
                                            Pw(:,:,ispin),ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin), &
                                            Fx(:,:,ispin),FxHF(:,:,ispin),Cx_choice,doNcentered)
     end do
@@ -338,7 +338,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
 !   Exchange energy
 
     do ispin=1,nspin
-      call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
+      call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, &
                                         Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin)&
                                         ,Cx_choice,doNcentered)
     end do
@@ -403,7 +403,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
 ! Compute individual energies from ensemble energy
 !------------------------------------------------------------------------
 
-  call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
+  call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,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,doNcentered)
 
diff --git a/src/eDFT/read_options_dft.f90 b/src/eDFT/read_options_dft.f90
index 72c7d44..86488f9 100644
--- a/src/eDFT/read_options_dft.f90
+++ b/src/eDFT/read_options_dft.f90
@@ -1,4 +1,4 @@
-subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, & 
+subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,nCC,aCC, & 
                             doNcentered,occnum,Cx_choice)
 
 ! Read DFT options
@@ -8,16 +8,17 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
   include 'parameters.h'
 
 ! Input variables
-  integer,intent(in)           :: nBas
+  integer,intent(in)            :: nBas
 
 ! Local variables
 
   integer                       :: iBas
   integer                       :: iEns
-  integer                       :: iParam
+  integer                       :: iCC
   character(len=1)              :: answer
   double precision,allocatable  :: nEl(:)
-  character(len=12)             :: x_func,c_func
+  character(len=12)             :: x_func
+  character(len=12)             :: c_func
 
 ! Output variables
 
@@ -28,8 +29,8 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
   integer,intent(out)           :: nEns
   logical,intent(out)           :: doNcentered
   double precision,intent(out)  :: wEns(maxEns)
-  double precision,intent(out)  :: aCC_w1(3)
-  double precision,intent(out)  :: aCC_w2(3)
+  integer,intent(out)           :: nCC
+  double precision,intent(out)  :: aCC(maxCC,maxEns-1)
   double precision,intent(out)  :: occnum(nBas,nspin,maxEns)
 
   integer,intent(out)           :: Cx_choice
@@ -394,24 +395,26 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
   
 ! Read parameters for weight-dependent functional
   read(1,*)
-  read(1,*) (aCC_w1(iParam),iParam=1,3)
-  read(1,*) (aCC_w2(iParam),iParam=1,3)
+  read(1,*) nCC
+  do iEns=2,nEns
+    read(1,*) (aCC(iCC,iEns-1),iCC=1,nCC)
+  end do
+
 ! Read choice of exchange coefficient
   read(1,*)
   read(1,*) Cx_choice
 
   write(*,*)'----------------------------------------------------------'
-  write(*,*)' parameters for w1-dependent exchange functional coefficient '
+  write(*,*)' Parameters for weight-dependent exchange functional      '
   write(*,*)'----------------------------------------------------------'
-  call matout(3,1,aCC_w1)
+  do iEns=2,nEns
+    write(*,'(A6,I2,A2)') 'State ',iEns,':'
+    do iCC=1,nCC
+      write(*,'(I2,F10.6)') iCC,aCC(iCC,iEns-1)
+    end do
+  end do
   write(*,*)
 
-  write(*,*)'----------------------------------------------------------'
-  write(*,*)' parameters for w2-dependent exchange functional coefficient '
-  write(*,*)'----------------------------------------------------------'
-  call matout(3,1,aCC_w2)
-  write(*,*) 
-
 ! Close file with options
 
   close(unit=1)
diff --git a/src/eDFT/unrestricted_exchange_energy.f90 b/src/eDFT/unrestricted_exchange_energy.f90
index f4db383..9b0f348 100644
--- a/src/eDFT/unrestricted_exchange_energy.f90
+++ b/src/eDFT/unrestricted_exchange_energy.f90
@@ -1,4 +1,4 @@
-subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF, &
+subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, &
                                         rho,drho,Ex,Cx_choice,doNcentered)
 
 ! Compute the exchange energy
@@ -13,8 +13,8 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,a
   logical,intent(in)            :: LDA_centered
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: wEns(nEns)
-  double precision,intent(in)   :: aCC_w1(3)
-  double precision,intent(in)   :: aCC_w2(3)
+  integer,intent(in)            :: nCC
+  double precision,intent(in)   :: aCC(nCC,nEns-1)
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   integer,intent(in)            :: nBas
@@ -43,7 +43,7 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,a
 
     case(1) 
 
-      call unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,&
+      call unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,&
                                             rho,Ex,Cx_choice,doNcentered)
 
 !   GGA functionals
@@ -62,7 +62,7 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,a
 
     case(4) 
 
-      call unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF, &
+      call unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, &
                                         rho,drho,Ex,Cx_choice)
 
   end select
diff --git a/src/eDFT/unrestricted_exchange_potential.f90 b/src/eDFT/unrestricted_exchange_potential.f90
index a52bca3..5c674e3 100644
--- a/src/eDFT/unrestricted_exchange_potential.f90
+++ b/src/eDFT/unrestricted_exchange_potential.f90
@@ -1,4 +1,4 @@
-subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P, & 
+subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, & 
                                            ERI,AO,dAO,rho,drho,Fx,FxHF,Cx_choice,doNcentered)
 
 ! Compute the exchange potential 
@@ -13,8 +13,8 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w
   logical,intent(in)            :: LDA_centered
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: wEns(nEns)
-  double precision,intent(in)   :: aCC_w1(3)
-  double precision,intent(in)   :: aCC_w2(3)
+  integer,intent(in)            :: nCC
+  double precision,intent(in)   :: aCC(nCC,nEns-1)
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   integer,intent(in)            :: nBas
@@ -50,7 +50,7 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w
 
     case(1) 
 
-      call unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,&
+      call unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Fx,&
                                                Cx_choice,doNcentered)
 
 !   GGA functionals
@@ -69,7 +69,7 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w
 
     case(4) 
 
-      call unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P, &
+      call unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, &
                                                   ERI,AO,dAO,rho,drho,Fx,FxHF,Cx_choice)
 
   end select
diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90
index 5acd636..247e94b 100644
--- a/src/eDFT/unrestricted_individual_energy.f90
+++ b/src/eDFT/unrestricted_individual_energy.f90
@@ -150,21 +150,6 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
 
   end do
 
-!------------------------------------------------------------------------
-! Checking Hartree contributions for each individual states
-!------------------------------------------------------------------------
-
-! print*,'Hartree contributions for each individual states'
-! print*,''
-! print*,''
-! print*,'EJ(aa,1)=',EJ(1,1),'EJ(ab,1)=',EJ(2,1),'EJ(bb,1)=',EJ(3,1)
-! print*,''
-! print*,'EJ(aa,2)=',EJ(1,2),'EJ(ab,2)=',EJ(2,2),'EJ(bb,2)=',EJ(3,2)
-! print*,''
-! print*,'EJ(aa,3)=',EJ(1,3),'EJ(ab,3)=',EJ(2,3),'EJ(bb,3)=',EJ(3,3)
-! print*,''
-
-
 !------------------------------------------------------------------------
 ! Individual exchange energy
 !------------------------------------------------------------------------
@@ -175,37 +160,10 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
                                                    Pw(:,:,ispin),P(:,:,ispin,iEns),rhow(:,ispin),drhow(:,:,ispin),           &  
                                                    rho(:,ispin,iEns),drho(:,:,ispin,iEns),Cx_choice,doNcentered,kappa(iEns), & 
                                                    Ex(ispin,iEns))
+
     end do
   end do
 
-!------------------------------------------------------------------------
-! Checking exchange contributions for each individual states
-!------------------------------------------------------------------------
-! print*,''
-! print*,''
-! print*,'Exchange contributions for each individual states'
-! print*,''
-! print*,''
-! print*,'Ex(aa,1) =' ,Ex(1,1),'Ex(bb,1) =' ,Ex(2,1)
-! print*,''
-! print*,'Ex(aa,2) =' ,Ex(1,2),'Ex(bb,2) =' ,Ex(2,2)
-! print*,''
-! print*,'Ex(aa,3) =' ,Ex(1,3),'Ex(bb,3) =' ,Ex(2,3)
-
-!------------------------------------------------------------------------
-! Checking number of alpha and beta electrons for each individual states
-!------------------------------------------------------------------------
-! print*,''
-! print*,''
-! print*,'Checking number of alpha and beta electrons for each individual states'
-! print*,''
-! print*,''
-! print*,'nEl(a,1) = ',electron_number(nGrid,weight,rho(:,1,1)),'nEl(b,1) = ',electron_number(nGrid,weight,rho(:,2,1))
-! print*,''
-! print*,'nEl(a,2) = ',electron_number(nGrid,weight,rho(:,1,2)),'nEl(b,2) = ',electron_number(nGrid,weight,rho(:,2,2))
-! print*,''
-! print*,'nEl(a,3) = ',electron_number(nGrid,weight,rho(:,1,3)),'nEl(b,3) = ',electron_number(nGrid,weight,rho(:,2,3))
-
 !------------------------------------------------------------------------
 ! Individual correlation energy
 !------------------------------------------------------------------------
diff --git a/src/eDFT/unrestricted_lda_exchange_energy.f90 b/src/eDFT/unrestricted_lda_exchange_energy.f90
index 4eb014d..4afa730 100644
--- a/src/eDFT/unrestricted_lda_exchange_energy.f90
+++ b/src/eDFT/unrestricted_lda_exchange_energy.f90
@@ -1,4 +1,4 @@
-subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,Cx_choice,doNcentered)
+subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rho,Ex,Cx_choice,doNcentered)
 
 ! Select LDA exchange functional
 
@@ -11,8 +11,8 @@ subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aC
   logical,intent(in)            :: LDA_centered
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: wEns(nEns)
-  double precision,intent(in)   :: aCC_w1(3)
-  double precision,intent(in)   :: aCC_w2(3)
+  integer,intent(in)            :: nCC
+  double precision,intent(in)   :: aCC(nCC,nEns-1)
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   double precision,intent(in)   :: rho(nGrid)
@@ -33,7 +33,7 @@ subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aC
 
     case (2)
 
-      call UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Cx_choice,doNcentered,Ex)
+      call UCC_lda_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex)
 
     case default
 
diff --git a/src/eDFT/unrestricted_lda_exchange_potential.f90 b/src/eDFT/unrestricted_lda_exchange_potential.f90
index ee7912c..21833fb 100644
--- a/src/eDFT/unrestricted_lda_exchange_potential.f90
+++ b/src/eDFT/unrestricted_lda_exchange_potential.f90
@@ -1,4 +1,4 @@
-subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx &
+subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Fx &
                                                ,Cx_choice,doNcentered)
 
 ! Select LDA correlation potential
@@ -13,8 +13,8 @@ subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1
   integer,intent(in)            :: DFA
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: wEns(nEns)
-  double precision,intent(in)   :: aCC_w1(3)
-  double precision,intent(in)   :: aCC_w2(3)
+  integer,intent(in)            :: nCC
+  double precision,intent(in)   :: aCC(nCC,nEns-1)
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   integer,intent(in)            :: nBas
@@ -37,7 +37,7 @@ subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1
 
     case (2)
 
-      call UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice,doNcentered)
+      call UCC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Fx,Cx_choice,doNcentered)
 
     case default