diff --git a/.gitignore b/.gitignore
index b2120b38..226118b4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,2 @@
 compile_commands.json
 doc/cpp2rst_generated
-build/
diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py
index f4aa22ed..14d16bf0 100644
--- a/python/triqs_dft_tools/sumk_dft.py
+++ b/python/triqs_dft_tools/sumk_dft.py
@@ -38,10 +38,11 @@ from warnings import warn
 from scipy import compress
 from scipy.optimize import minimize, newton, brenth
 
-def compute_DC_from_density(N_tot, U, J, N_spin = None,  n_orbitals=5,  method='cFLL'):
+def compute_DC_from_density(N_tot, U, J, N_spin = None,  n_orbitals=5,  method='sFLL'):
     """
-    Computes the double counting correction in terms using various methods
-
+    Computes the double counting correction using various methods.
+    For FLL and AMF DC the notations and equations from doi.org/10.1038/s41598-018-27731-4
+    are used, whereas for the Held DC the definitions from doi.org/10.1080/00018730701619647 are used.
 
     Parameters
     ----------
@@ -69,45 +70,48 @@ def compute_DC_from_density(N_tot, U, J, N_spin = None,  n_orbitals=5,  method='
     List of floats:
         -   DC_val: double counting potential
         -   E_val: double counting energy
+
+
+    todo: See whether to move this to TRIQS directly instead of dft_tools
     """
     if N_spin is not None:
         N_spin2 = N_tot-N_spin
         Mag = N_spin - N_spin2
     L_orbit = (n_orbitals-1)/2 
 
-    match method:
-        case 'cFLL':
-            E_val = 0.5 * U * N_tot * (N_tot-1) - 0.5 * J * N_tot * (N_tot*0.5-1)
-            DC_val = U * (N_tot-0.5) - J *(N_tot*0.5-0.5)
+    if method == 'cFLL':
+        E_val = 0.5 * U * N_tot * (N_tot-1) - 0.5 * J * N_tot * (N_tot*0.5-1)
+        DC_val = U * (N_tot-0.5) - J *(N_tot*0.5-0.5)
 
-        case 'sFLL':
-            assert N_spin is not None, "Spin density not given"
-            E_val = 0.5 * U * N_tot * (N_tot-1) - 0.5 * J * N_tot * (N_tot*0.5-1) - 0.25 * J * Mag**2
-            DC_val = U * (N_tot-0.5) - J *(N_spin-0.5)
+    elif method == 'sFLL':
+        assert N_spin is not None, "Spin density not given"
+        E_val = 0.5 * U * N_tot * (N_tot-1) - 0.5 * J * N_tot * (N_tot*0.5-1) - 0.25 * J * Mag**2
+        DC_val = U * (N_tot-0.5) - J *(N_spin-0.5)
         
-        case 'cAMF':
-            E_val = +0.5 * U * N_tot **2
-            E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot**2
-            DC_val = U * N_tot - 0.5*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot
+    elif method == 'cAMF':
+        E_val = +0.5 * U * N_tot **2
+        E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot**2
+        DC_val = U * N_tot - 0.5*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot
             
 
-        case 'sAMF':
-            assert N_spin is not None, "Spin density not given"
-            E_val = 0.5 * U * N_tot **2
-            E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot**2
-            E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*Mag**2
-            DC_val = U * N_tot  - (U+2*L_orbit*J)/(2*L_orbit+1)*N_spin
+    elif method == 'sAMF':
+        assert N_spin is not None, "Spin density not given"
+        E_val = 0.5 * U * N_tot **2
+        E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot**2
+        E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*Mag**2
+        DC_val = U * N_tot  - (U+2*L_orbit*J)/(2*L_orbit+1)*N_spin
         
-        case 'cHeld':
-            U_mean = (U + (n_orbitals-1)*(U-2*J)+(n_orbitals-1)*(U-3*J))/(2*n_orbitals-1)
-            E_val = 0.5 * U_mean * N_tot * (N_tot - 1)
-            DC_val = U_mean * (N_tot-0.5)
+    elif method == 'cHeld':
+        # Valid for a Kanamori-type Hamiltonian where U'=U-2J
+        U_mean = (U + (n_orbitals-1)*(U-2*J)+(n_orbitals-1)*(U-3*J))/(2*n_orbitals-1)
+        E_val = 0.5 * U_mean * N_tot * (N_tot - 1)
+        DC_val = U_mean * (N_tot-0.5)
         
-        case 'sHeld':
-            raise ValueError(f"Method sHeld not yet implemented")
+    elif method == 'sHeld':
+        raise ValueError(f"Method sHeld not yet implemented")
         
-        case _:
-            raise ValueError(f"DC type {method} not supported")
+    else:
+        raise ValueError(f"DC type {method} not supported")
 
     mpi.report(f"DC potential computed using the {method} method, V_DC = {DC_val:.6f} eV")
     mpi.report(f"E_DC using the {method} method, E_DC = {E_val:.6f} eV")
@@ -1735,7 +1739,7 @@ class SumkDFT(object):
         self.dc_energ = dc_energ
 
     def calc_dc(self, dens_mat, orb=0, U_interact=None, J_hund=None,
-                use_dc_formula=0, use_dc_value=None, transform=True, legacy=True):
+                use_dc_formula=0, use_dc_value=None, transform=True):
         r"""
         Calculate and set the double counting corrections.
 
@@ -1767,8 +1771,13 @@ class SumkDFT(object):
             Value of interaction parameter `U`.
         J_hund : float, optional
             Value of interaction parameter `J`.
-        use_dc_formula : int, optional
-            Type of double-counting correction (see description).
+        use_dc_formula : int or string, optional 
+            Type of double-counting correction (see description of `compute_DC_from_density` above).
+            There is an interface with the legacy implementation which allows for the old convention:
+            * 0 -> 'sFLL'  spin dependent fully localized limit
+            * 1 -> 'cHeld' spin independent Held formula
+            * 2 -> 'sAMF'  spin dependent around-mean field approximation
+
         use_dc_value : float, optional
             Value of the double-counting correction. If specified
             `U_interact`, `J_hund` and `use_dc_formula` are ignored.
@@ -1807,68 +1816,21 @@ class SumkDFT(object):
                 dim //= 2
             
             if use_dc_value is None:
-                if legacy:
-
-                    if U_interact is None and J_hund is None:
-                        raise ValueError("set_dc: either provide U_interact and J_hund or set use_dc_value to dc value.")
-
-                    if use_dc_formula == 0:  # FLL
-
-                        self.dc_energ[icrsh] = U_interact / \
-                            2.0 * Ncrtot * (Ncrtot - 1.0)
-                        for sp in spn:
-                            Uav = U_interact * (Ncrtot - 0.5) - \
-                                J_hund * (Ncr[sp] - 0.5)
-                            self.dc_imp[icrsh][sp] *= Uav
-                            self.dc_energ[icrsh] -= J_hund / \
-                                len(spn) * (Ncr[sp]) * (Ncr[sp] - 1.0)
-                            mpi.report(
-                                "DC for shell %(icrsh)i and block %(sp)s = %(Uav)f" % locals())
-
-                    elif use_dc_formula == 1:  # Held's formula, with U_interact the interorbital onsite interaction
-
-                        self.dc_energ[icrsh] = (U_interact + (dim - 1) * (U_interact - 2.0 * J_hund) + (
-                            dim - 1) * (U_interact - 3.0 * J_hund)) / (2 * dim - 1) / 2.0 * Ncrtot * (Ncrtot - 1.0)
-                        for sp in spn:
-                            Uav = (U_interact + (dim - 1) * (U_interact - 2.0 * J_hund) + (dim - 1)
-                                   * (U_interact - 3.0 * J_hund)) / (2 * dim - 1) * (Ncrtot - 0.5)
-                            self.dc_imp[icrsh][sp] *= Uav
-                            mpi.report(
-                                "DC for shell %(icrsh)i and block %(sp)s = %(Uav)f" % locals())
-
-                    elif use_dc_formula == 2:  # AMF
-
-                        self.dc_energ[icrsh] = 0.5 * U_interact * Ncrtot * Ncrtot
-                        for sp in spn:
-                            Uav = U_interact * \
-                                (Ncrtot - Ncr[sp] / dim) - \
-                                J_hund * (Ncr[sp] - Ncr[sp] / dim)
-                            self.dc_imp[icrsh][sp] *= Uav
-                            self.dc_energ[
-                                icrsh] -= (U_interact + (dim - 1) * J_hund) / dim / len(spn) * Ncr[sp] * Ncr[sp]
-                            mpi.report(
-                                "DC for shell %(icrsh)i and block %(sp)s = %(Uav)f" % locals())
-
-                    mpi.report("DC energy for shell %s = %s" %
-                               (icrsh, self.dc_energ[icrsh]))
+                #For legacy compatibility 
+                if use_dc_formula == 0:
+                    mpi.report(f"Detected {use_dc_formula=}, changing to sFLL")
+                    use_dc_formula = "sFLL"
+                if use_dc_formula == 1:
+                    mpi.report(f"Detected {use_dc_formula=}, changing to cHeld")
+                    use_dc_formula = "cHeld"
+                if use_dc_formula == 2:
+                    mpi.report(f"Detected {use_dc_formula=}, changing to sAMF")
+                    use_dc_formula = "sAMF"
                 
-                else:
-                    #For legacy compatibility 
-                    match use_dc_formula:
-                        case 0:
-                            mpi.report(f"Detected {use_dc_formula=}, changing to sFLL")
-                            use_dc_formula = "sFLL"
-                        case 1:
-                            mpi.report(f"Detected {use_dc_formula=}, changing to cHeld")
-                            use_dc_formula = "cHeld"
-                        case 2:
-                            mpi.report(f"Detected {use_dc_formula=}, changing to sAMF")
-                            use_dc_formula = "sAMF"
-                    
-                    for sp in spn:
-                        DC_val, E_val = compute_DC_from_density(N_tot=Ncrtot,U=U_interact, J=J_hund, n_orbitals=dim, N_spin=Ncr[sp], method=use_dc_formula)
-                        self.dc_imp[icrsh][sp] *= DC_val
-                        self.dc_energ[icrsh] -= E_val
+                for sp in spn:
+                    DC_val, E_val = compute_DC_from_density(N_tot=Ncrtot,U=U_interact, J=J_hund, n_orbitals=dim, N_spin=Ncr[sp], method=use_dc_formula)
+                    self.dc_imp[icrsh][sp] *= DC_val
+                    self.dc_energ[icrsh] = E_val
 
 
 
@@ -2048,7 +2010,6 @@ class SumkDFT(object):
     def calc_mu(self, precision=0.01, broadening=None, delta=0.5, max_loops=100, method="dichotomy"):
         r"""
         Searches for the chemical potential that gives the DFT total charge.
-        A simple bisection method is used.
 
         Parameters
         ----------
@@ -2061,7 +2022,7 @@ class SumkDFT(object):
         max_loops : int, optional
                     Number of dichotomy loops maximally performed.
         
-        max_loops : string, optional
+        method : string, optional
                     Type of optimization used:
                         * dichotomy: usual bisection algorithm from the TRIQS library
                         * newton: newton method, faster convergence but more unstable 
@@ -2116,52 +2077,50 @@ class SumkDFT(object):
 
             mpi.report("Trying out mu = {}".format(str(mu)))
             calc_dens = self.total_density(mu=mu, broadening=broadening).real - density
-            mpi.report("Delta to target density = {}".format(str(calc_dens)))
+            mpi.report(f"Target density = {density}; Delta to target = {calc_dens}")
             return calc_dens
         
         #check for lowercase matching for the method variable
-        match method.lower():
+        if method.lower()=="dichotomy":
+            mpi.report("\nSUMK calc_mu: Using dichtomy adjustment to find chemical potential\n")
+            self.chemical_potential = dichotomy.dichotomy(function=F_bisection,
+                                                          x_init=self.chemical_potential, y_value=density,
+                                                          precision_on_y=precision, delta_x=delta, max_loops=max_loops,
+                                                          x_name="Chemical Potential", y_name="Total Density",
+                                                          verbosity=3)[0]
+        elif method.lower()=="newton":
+            mpi.report("\nSUMK calc_mu: Using Newton method to find chemical potential\n")
+            self.chemical_potential =                newton(func=F_optimize,
+                                                          x0=self.chemical_potential,
+                                                          tol=precision, maxiter=max_loops,
+                                                          )
 
-            case "dichotomy":
-                mpi.report("\nSUMK calc_mu: Using dichtomy adjustment to find chemical potential\n")
-                self.chemical_potential = dichotomy.dichotomy(function=F_bisection,
-                                                              x_init=self.chemical_potential, y_value=density,
-                                                              precision_on_y=precision, delta_x=delta, max_loops=max_loops,
-                                                              x_name="Chemical Potential", y_name="Total Density",
-                                                              verbosity=3)[0]
-            case 'newton':
-                mpi.report("\nSUMK calc_mu: Using Newton method to find chemical potential\n")
-                self.chemical_potential =                newton(func=F_optimize,
-                                                              x0=self.chemical_potential,
-                                                              tol=precision, maxiter=max_loops,
-                                                              )
-
-            case 'brent':
-                mpi.report("\nSUMK calc_mu: Using Brent method to find chemical potential")
-                mpi.report("SUMK calc_mu: Finding bounds \n")
-                
-                mu_guess_0, mu_guess_1  =  find_bounds(function=F_optimize,
-                                                      x_init=self.chemical_potential,
-                                                      delta_x=delta, max_loops=max_loops,
-                                                      )
-                mu_guess_1 += 0.01 #scrambles higher lying interval to avoid getting stuck
-                mpi.report("\nSUMK calc_mu: Searching root with Brent method\n")
-                self.chemical_potential = brenth(f=F_optimize,
-                                                   a=mu_guess_0,
-                                                   b=mu_guess_1,
-                                                   xtol=precision, maxiter=max_loops,
-                                                   )
+        elif method.lower()=="brent":
+            mpi.report("\nSUMK calc_mu: Using Brent method to find chemical potential")
+            mpi.report("SUMK calc_mu: Finding bounds \n")
             
-            case _:
-                raise ValueError(
-                        f"SUMK calc_mu: The method selected: {method}, is not implemented\n",
-                        """
-                        Please check for typos or select one of the following:
-                            * dichotomy: usual bisection algorithm from the TRIQS library
-                            * newton: newton method, fastest convergence but more unstable 
-                            * brent: finds bounds and proceeds with hyperbolic brent method, a compromise between speed and ensuring convergence
-                        """
-                        )
+            mu_guess_0, mu_guess_1  =  find_bounds(function=F_optimize,
+                                                  x_init=self.chemical_potential,
+                                                  delta_x=delta, max_loops=max_loops,
+                                                  )
+            mu_guess_1 += 0.01 #scrambles higher lying interval to avoid getting stuck
+            mpi.report("\nSUMK calc_mu: Searching root with Brent method\n")
+            self.chemical_potential = brenth(f=F_optimize,
+                                               a=mu_guess_0,
+                                               b=mu_guess_1,
+                                               xtol=precision, maxiter=max_loops,
+                                               )
+            
+        else:
+            raise ValueError(
+                    f"SUMK calc_mu: The selected method: {method}, is not implemented\n",
+                    """
+                    Please check for typos or select one of the following:
+                        * dichotomy: usual bisection algorithm from the TRIQS library
+                        * newton: newton method, fastest convergence but more unstable 
+                        * brent: finds bounds and proceeds with hyperbolic brent method, a compromise between speed and ensuring convergence
+                    """
+                    )
 
         return self.chemical_potential
 
@@ -2414,7 +2373,8 @@ class SumkDFT(object):
                             f.write("\n")
 
         elif dm_type == 'qe':
-            assert self.SP == 0, "Spin-polarized density matrix is not implemented"
+            if self.SP == 0:
+                mpi.report("SUMK calc_density_correction: WARNING! Averaging out spin-polarized correction in the density channel")
 
             subgrp = 'dft_update'
             delta_N = np.zeros([self.n_k, max(self.n_orbitals[:,0]), max(self.n_orbitals[:,0])], dtype=complex)
diff --git a/test/python/NiO.ref.h5 b/test/python/NiO.ref.h5
new file mode 100755
index 00000000..13140104
Binary files /dev/null and b/test/python/NiO.ref.h5 differ
diff --git a/test/python/NiO.ref.h5_copy b/test/python/NiO.ref.h5_copy
new file mode 100755
index 00000000..24d918ca
Binary files /dev/null and b/test/python/NiO.ref.h5_copy differ
diff --git a/test/python/SrVO3.ref.h5 b/test/python/SrVO3.ref.h5
index 4150e050..70f848b0 100644
Binary files a/test/python/SrVO3.ref.h5 and b/test/python/SrVO3.ref.h5 differ
diff --git a/test/python/SrVO3.ref.h5_copy b/test/python/SrVO3.ref.h5_copy
new file mode 100644
index 00000000..c2a24b63
Binary files /dev/null and b/test/python/SrVO3.ref.h5_copy differ
diff --git a/test/python/dc_test.py b/test/python/dc_test.py
new file mode 100755
index 00000000..ba003355
--- /dev/null
+++ b/test/python/dc_test.py
@@ -0,0 +1,168 @@
+##########################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2023 by A. Carta
+#
+# TRIQS is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License along with
+# TRIQS. If not, see <http://www.gnu.org/licenses/>.
+#
+##########################################################################
+
+
+
+import numpy as np
+from triqs_dft_tools.sumk_dft import *
+from triqs.gf import  *
+from h5 import HDFArchive
+from triqs.operators.util import *
+import triqs.utility.mpi as mpi
+
+
+Uval = 5
+Jval = 0.3
+
+method_dict = {
+        "FLL":{
+            "numbering_convention":0,
+            "new_convention":"cFLL"
+            },
+        "AMF":{
+            "numbering_convention":2,
+            "new_convention":"cAMF"
+            },
+        "Held":{
+            "numbering_convention":1,
+            "new_convention":"cHeld"
+            },
+        }
+
+
+
+
+def test_dc(method, method_dict, dens, Uval, Jval, filename):
+
+    dc_no = method_dict[method]["numbering_convention"]
+    dc_string = method_dict[method]["new_convention"]
+    
+    mpi.report("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX")
+    mpi.report(f"\n Testing interface {method} \n")
+    mpi.report("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX")
+
+    mpi.report("\nTesting legacy compatibility layer:\n")
+    SK_compat.calc_dc(dens_mat=dens[0], U_interact=Uval, J_hund=Jval, use_dc_formula=dc_no )
+    mpi.report("Up DC matrix:")
+    mpi.report(SK_compat.dc_imp[0]['up'])
+    mpi.report(f"Double counting energy = {SK_compat.dc_energ} ")
+
+    mpi.report("\nTesting new dc interface:\n")
+    SK_new.calc_dc(dens_mat=dens[0], U_interact=Uval, J_hund=Jval, use_dc_formula=dc_string)
+    mpi.report("Up DC matrix:")
+    mpi.report(SK_new.dc_imp[0]['up'])
+    mpi.report(f"Double counting energy = {SK_new.dc_energ} ")
+    
+    # Load previously computed DC values from h5 archive 
+    R = HDFArchive(f'./{filename}', 'r')
+    dc_comp = R[f'DC_{method}_benchmark']['dc_imp']
+    en_comp = R[f'DC_{method}_benchmark']['dc_energ']
+    del R
+
+    mpi.report(f"\nAsserting comparison for method: {method}")
+    assert np.allclose(SK_compat.dc_imp[0]['up'], dc_comp, atol=1e-12), f"Assertion failed comparing legacy Vdc to reference, method: {method} "
+    assert np.allclose(SK_compat.dc_energ, en_comp, atol=1e-12), f"Assertion failed comparing legacy energy to reference, method {method} "
+    assert np.allclose(SK_new.dc_imp[0]['up'], dc_comp, atol=1e-12), f"Assertion failed comparing Vdc to reference, method: {method} "
+    assert np.allclose(SK_new.dc_energ, en_comp, atol=1e-12), f"Assertion failed comparing energy to reference, method: {method} "
+    mpi.report("Comparison with stored DC values successfull!\n")
+    
+
+
+
+# %% 5 orbitals testing
+
+mpi.report("\n############################################")
+mpi.report("############################################")
+mpi.report(f"\n \n Starting tests for 5 orbitals \n \n")
+mpi.report("############################################")
+mpi.report("############################################")
+dft_filename = "./NiO.ref"
+use_blocks = False
+SK_compat = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks)
+SK_compat.set_mu(13.9)
+SK_new    = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks)
+SK_new.set_mu(13.9)
+
+
+icrsh = 0
+dens = SK_compat.density_matrix()
+
+with np.printoptions(precision=5):
+    for key in dens[0].keys():
+        mpi.report(f"{key} channel")
+        mpi.report(dens[0][key].real)
+
+N_up = np.trace(dens[0]['up'].real)
+N_down = np.trace(dens[0]['down'].real)
+N_tot = N_up + N_down
+
+mpi.report(f"{N_up=} ,{N_down=}, {N_tot=}\n")
+
+for method in ["FLL", "AMF", "Held"]: 
+    test_dc(method, method_dict, dens, Uval, Jval, filename = f"{dft_filename}.h5")
+
+    #in case implementation changes, to write new testing data into archive
+    #R = HDFArchive('./NiO.ref.h5', 'a')
+    #R.create_group(f'DC_{method}_benchmark')
+    #R[f'DC_{method}_benchmark']['dc_imp']= SK_new.dc_imp[0]['up']
+    #R[f'DC_{method}_benchmark']['dc_energ']= SK_new.dc_energ
+    #del R
+
+
+# 3 orbital testing
+
+mpi.report("############################################")
+mpi.report("############################################")
+mpi.report(f"\n \n Starting tests for 3 orbitals \n \n")
+mpi.report("############################################")
+mpi.report("############################################")
+dft_filename = "./SrVO3.ref"
+use_blocks = False
+SK_compat = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks)
+SK_new    = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks)
+
+icrsh = 0
+dens = SK_compat.density_matrix()
+
+
+with np.printoptions(precision=5):
+    for key in dens[0].keys():
+        mpi.report(f"{key} channel")
+        mpi.report(dens[0][key].real)
+
+N_up = np.trace(dens[0]['up'].real)
+N_down = np.trace(dens[0]['down'].real)
+N_tot = N_up + N_down
+
+mpi.report(f"{N_up=} ,{N_down=}, {N_tot=}\n")
+
+Uval = 5
+Jval = 0.3
+
+for method in ["FLL", "AMF", "Held"]: 
+    test_dc(method, method_dict, dens, Uval, Jval, filename = f"{dft_filename}.h5" )
+    
+    #in case implementation changes, to write new testing data into archive
+    #R = HDFArchive(f'./{dft_filename}.h5', 'a')
+    #R.create_group(f'DC_{method}_benchmark')
+    #R[f'DC_{method}_benchmark']['dc_imp']= SK_new.dc_imp[0]['up']
+    #R[f'DC_{method}_benchmark']['dc_energ']= SK_new.dc_energ
+    #del R