[feat] improved standard behavior of block struct (#248)

* previously the default gf_struct_solver had keys up / down,
inconsistent with the default behavior after analyse_block_structure was
run: up_0 / down_0. Now the default solver structure always has the _0
in the key.
* old behavior resulted in error when analyse_block_structure was called
twice
* fixed analyse block structure tests with new changes
* to correctly use analyse_block_structure use now
extract_G_loc(transform_to_solver_blocks=False)
* changed density_matrix function to use directly extract_G_loc() if
using_gf is selected.
* print deprecation warning in density_matrix, same can be achieved via
extract_G_loc and [G.density() for G in Gloc]
* new function density_matrix_using_point_integration()
* enforce in analyse block structure that input dm or G is list with
length of n_corr_shells
* correct doc string for how include_shells are given
* fixed other tests accordingly
* fixed small bug in initial block structure regarding length of lists

python/triqs_dft_tools/sumk_dft.py

@@ -168,8 +168,8 @@ class SumkDFT(object):
             # blocks possible
             self.gf_struct_sumk = [[(sp, self.corr_shells[icrsh]['dim']) for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]]
                                    for icrsh in range(self.n_corr_shells)]
-            # First set a standard gf_struct solver:
+            # First set a standard gf_struct solver (add _0 here for consistency with analyse_block_structure):
-            self.gf_struct_solver = [dict([(sp, self.corr_shells[self.inequiv_to_corr[ish]]['dim'])
+            self.gf_struct_solver = [dict([(sp+'_0', self.corr_shells[self.inequiv_to_corr[ish]]['dim'])
                                            for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]])
                                      for ish in range(self.n_inequiv_shells)]
             # Set standard (identity) maps from gf_struct_sumk <->
@@ -180,12 +180,12 @@ class SumkDFT(object):
                                          for ish in range(self.n_inequiv_shells)]
             for ish in range(self.n_inequiv_shells):
                 for block, inner_dim in self.gf_struct_sumk[self.inequiv_to_corr[ish]]:
-                    self.solver_to_sumk_block[ish][block] = block
+                    self.solver_to_sumk_block[ish][block+'_0'] = block
                     for inner in range(inner_dim):
                         self.sumk_to_solver[ish][
-                            (block, inner)] = (block, inner)
+                            (block, inner)] = (block+'_0', inner)
                         self.solver_to_sumk[ish][
-                            (block, inner)] = (block, inner)
+                            (block+'_0', inner)] = (block, inner)
             # assume no shells are degenerate
             self.deg_shells = [[] for ish in range(self.n_inequiv_shells)]
 
@@ -862,8 +862,8 @@ class SumkDFT(object):
                     If the difference between density matrix / hloc elements is below threshold,
                     they are considered to be equal.
         include_shells : list of integers, optional
-                         List of correlated shells to be analysed.
+                         List of inequivalent shells to be analysed.
-                         If include_shells is not provided all correlated shells will be analysed.
+                         If include_shells is not provided all inequivalent shells will be analysed.
         dm : list of dict, optional
              List of density matrices from which block stuctures are to be analysed.
              Each density matrix is a dict {block names: 2d numpy arrays} for each correlated shell.
@@ -874,14 +874,11 @@ class SumkDFT(object):
                If not provided, it will be calculated using eff_atomic_levels.
         """
 
-        self.gf_struct_solver = [{} for ish in range(self.n_inequiv_shells)]
-        self.sumk_to_solver = [{} for ish in range(self.n_inequiv_shells)]
-        self.solver_to_sumk = [{} for ish in range(self.n_inequiv_shells)]
-        self.solver_to_sumk_block = [{}
-                                     for ish in range(self.n_inequiv_shells)]
-
         if dm is None:
-            dm = self.density_matrix(method='using_point_integration')
+            warn("WARNING: No density matrix given. Calculating density matrix with default parameters. This will be deprecated in future releases.")
+            dm = self.density_matrix(method='using_gf', transform_to_solver_blocks=False)
+
+        assert len(dm) == self.n_corr_shells, "The number of density matrices must be equal to the number of correlated shells."
         dens_mat = [dm[self.inequiv_to_corr[ish]]
                     for ish in range(self.n_inequiv_shells)]
         if hloc is None:
@@ -890,8 +887,13 @@ class SumkDFT(object):
         if include_shells is None:
             include_shells = list(range(self.n_inequiv_shells))
         for ish in include_shells:
+            self.gf_struct_solver[ish] = {}
+            self.sumk_to_solver[ish] = {}
+            self.solver_to_sumk[ish] = {}
+            self.solver_to_sumk_block[ish] = {}
 
             for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
+                assert sp in dens_mat[ish], f"The density matrix does not contain the block {sp}. Is the input dm given in sumk block structure?"
                 n_orb = self.corr_shells[self.inequiv_to_corr[ish]]['dim']
                 # gives an index list of entries larger that threshold
                 dmbool = (abs(dens_mat[ish][sp]) > threshold)
@@ -1049,8 +1051,8 @@ class SumkDFT(object):
                     If the difference between matrix elements is below threshold,
                     they are considered to be equal.
         include_shells : list of integers, optional
-                         List of correlated shells to be analysed.
+                         List of inequivalent shells to be analysed.
-                         If include_shells is not provided all correlated shells will be analysed.
+                         If include_shells is not provided all inequivalent shells will be analysed.
         analyse_deg_shells : bool
                              Whether to call the analyse_deg_shells function
                              after having finished the block structure analysis
@@ -1062,29 +1064,26 @@ class SumkDFT(object):
         """
 
         assert isinstance(G, list), "G must be a list (with elements for each correlated shell)"
+        assert len(G) == self.n_corr_shells, "G must have one element for each correlated shell, run extract_G_loc with transform_to_solver_blocks=False to get the correct G"
 
         gf = self._get_hermitian_quantity_from_gf(G)
 
-        # initialize the variables
-        self.gf_struct_solver = [{} for ish in range(self.n_inequiv_shells)]
-        self.sumk_to_solver = [{} for ish in range(self.n_inequiv_shells)]
-        self.solver_to_sumk = [{} for ish in range(self.n_inequiv_shells)]
-        self.solver_to_sumk_block = [{}
-                                     for ish in range(self.n_inequiv_shells)]
-
-        # the maximum value of each matrix element of each block and shell
-        max_gf = [{name:np.max(np.abs(g.data),0) for name, g in gf[ish]} for ish in range(self.n_inequiv_shells)]
-
         if include_shells is None:
             # include all shells
             include_shells = list(range(self.n_inequiv_shells))
 
         for ish in include_shells:
+            self.gf_struct_solver[ish] = {}
+            self.sumk_to_solver[ish] = {}
+            self.solver_to_sumk[ish] = {}
+            self.solver_to_sumk_block[ish] = {}
             for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
+                assert sp in gf[self.inequiv_to_corr[ish]].indices, f"The Green's function does not contain the block {sp}. Is the input G given in sumk block structure?"
                 n_orb = self.corr_shells[self.inequiv_to_corr[ish]]['dim']
 
                 # gives an index list of entries larger that threshold
-                maxgf_bool = (abs(max_gf[ish][sp]) > threshold)
+                max_gf = np.max(np.abs(gf[self.inequiv_to_corr[ish]][sp].data),0)
+                maxgf_bool = (max_gf > threshold)
 
                 # Determine off-diagonal entries in upper triangular part of the
                 # Green's function
@@ -1455,21 +1454,12 @@ class SumkDFT(object):
 
         return trans
 
-
+    def density_matrix_using_point_integration(self):
-    def density_matrix(self, method='using_gf'):
+        """
-        """Calculate density matrices in one of two ways.
+        Calculate density matrices using point integration: Only works for
-
+        diagonal hopping matrix (true in wien2k). Consider using extract_G_loc
-        Parameters
+        together with [G.density() for G in Gloc] instead. Returned density
-        ----------
+        matrix is always given in SumK block structure.
-        method : string, optional
-
-                 - if 'using_gf': First get lattice gf (g_loc is not set up), then density matrix.
-                                  It is useful for Hubbard I, and very quick.
-                                  No assumption on the hopping structure is made (ie diagonal or not).
-                 - if 'using_point_integration': Only works for diagonal hopping matrix (true in wien2k).
-
-        beta : float, optional
-               Inverse temperature.
 
         Returns
         -------
@@ -1483,19 +1473,9 @@ class SumkDFT(object):
                     [self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], complex)
 
         ikarray = np.array(list(range(self.n_k)))
-        for ik in mpi.slice_array(ikarray):
-
-            if method == "using_gf":
-
-                G_latt_iw = self.lattice_gf(ik=ik, mu=self.chemical_potential)
-                G_latt_iw *= self.bz_weights[ik]
-                dm = G_latt_iw.density()
-                MMat = [dm[sp] for sp in self.spin_block_names[self.SO]]
-
-            elif method == "using_point_integration":
-
         ntoi = self.spin_names_to_ind[self.SO]
         spn = self.spin_block_names[self.SO]
+        for ik in mpi.slice_array(ikarray):
             dims = {sp:self.n_orbitals[ik, ntoi[sp]] for sp in spn}
             MMat = [np.zeros([dims[sp], dims[sp]], complex) for sp in spn]
 
@@ -1509,9 +1489,6 @@ class SumkDFT(object):
                     else:
                         MMat[isp][inu, inu] = 0.0
 
-            else:
-                raise ValueError("density_matrix: the method '%s' is not supported." % method)
-
             for icrsh in range(self.n_corr_shells):
                 for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]):
                     ind = self.spin_names_to_ind[
@@ -1519,10 +1496,6 @@ class SumkDFT(object):
                     dim = self.corr_shells[icrsh]['dim']
                     n_orb = self.n_orbitals[ik, ind]
                     projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb]
-                    if method == "using_gf":
-                        dens_mat[icrsh][sp] += np.dot(np.dot(projmat, MMat[isp]),
-                                                         projmat.transpose().conjugate())
-                    elif method == "using_point_integration":
                     dens_mat[icrsh][sp] += self.bz_weights[ik] * np.dot(np.dot(projmat, MMat[isp]),
                                                                                projmat.transpose().conjugate())
 
@@ -1546,6 +1519,55 @@ class SumkDFT(object):
 
         return dens_mat
 
+
+    def density_matrix(self, method='using_gf', mu=None, with_Sigma=True, with_dc=True, broadening=None,
+                      transform_to_solver_blocks=True, show_warnings=True):
+        """Calculate density matrices in one of two ways.
+
+        Parameters
+        ----------
+        method : string, optional
+
+                 - if 'using_gf': First get lattice gf (g_loc is not set up), then density matrix.
+                                  It is useful for Hubbard I, and very quick.
+                                  No assumption on the hopping structure is made (ie diagonal or not).
+                 - if 'using_point_integration': Only works for diagonal hopping matrix (true in wien2k).
+        mu : real, optional
+            Input chemical potential. If not provided the value of self.chemical_potential is used as mu.
+        with_Sigma : boolean, optional
+            If True then the local GF is calculated with the self-energy self.Sigma_imp.
+        with_dc : boolean, optional
+            If True then the double-counting correction is subtracted from the self-energy in calculating the GF.
+        broadening : float, optional
+            Imaginary shift for the axis along which the real-axis GF is calculated.
+            If not provided, broadening will be set to double of the distance between mesh points in 'mesh'.
+            Only relevant for real-frequency GF.
+        transform_to_solver_blocks : bool, optional
+            If True (default), the returned G_loc will be transformed to the block structure ``gf_struct_solver``,
+            else it will be in ``gf_struct_sumk``.
+        show_warnings : bool, optional
+            Displays warning messages during transformation
+            (Only effective if transform_to_solver_blocks = True
+
+        Returns
+        -------
+        dens_mat : list of dicts
+                   Density matrix for each spin in each correlated shell.
+        """
+
+        if method == "using_gf":
+            warn("WARNING: density_matrix: method 'using_gf' is deprecated. Use 'extract_G_loc' instead.")
+            Gloc = self.extract_G_loc(mu, with_Sigma, with_dc, broadening,
+                                      transform_to_solver_blocks, show_warnings)
+            dens_mat = [G.density() for G in Gloc]
+        elif method == "using_point_integration":
+            warn("WARNING: density_matrix: method 'using_point_integration' is deprecated. Use 'density_matrix_using_point_integration' instead. All additionally provided arguments are ignored.")
+            dens_mat = self.density_matrix_using_point_integration()
+        else:
+           raise ValueError("density_matrix: the method '%s' is not supported." % method)
+
+        return dens_mat
+
     # For simple dft input, get crystal field splittings.
     def eff_atomic_levels(self):
         r"""

test/python/analyse_block_structure_from_gf.py

@@ -1,9 +1,9 @@
-from triqs.gf import *
+from triqs.gf import MeshImFreq, inverse, Omega
 from triqs_dft_tools.sumk_dft import SumkDFT
 from scipy.linalg import expm
 import numpy as np
 from triqs.utility.comparison_tests import assert_gfs_are_close, assert_arrays_are_close, assert_block_gfs_are_close
-from h5 import *
+from h5 import HDFArchive
 import itertools
 
 # The full test checks all different possible combinations of conjugated
@@ -22,8 +22,7 @@ mesh = MeshImFreq(40, 'Fermion', 1025)
 SK = SumkDFT(hdf_file='SrIrO3_rot.h5', mesh=mesh)
 Sigma = SK.block_structure.create_gf(mesh=mesh)
 SK.put_Sigma([Sigma])
-G = SK.extract_G_loc()
+G = SK.extract_G_loc(transform_to_solver_blocks=False)
-
 # the original block structure
 block_structure1 = SK.block_structure.copy()
 G_new = SK.analyse_block_structure_from_gf(G)
@@ -36,12 +35,11 @@ with HDFArchive('analyse_block_structure_from_gf.out.h5','w') as ar:
     ar['bs2'] = block_structure2
 
 # check whether the block structure is the same as in the reference
-with HDFArchive('analyse_block_structure_from_gf.out.h5','r') as ar,\
+with HDFArchive('analyse_block_structure_from_gf.out.h5', 'r') as ar, HDFArchive('analyse_block_structure_from_gf.ref.h5', 'r') as ar2:
-     HDFArchive('analyse_block_structure_from_gf.ref.h5','r') as ar2:
     assert ar['bs1'] == ar2['bs1'], 'bs1 not equal'
     a1 = ar['bs2']
     a2 = ar2['bs2']
-    assert a1==block_structure2, "writing/reading block structure incorrect"
+    assert a1 == block_structure2, 'writing/reading block structure incorrect'
     # we set the deg_shells to None because the transformation matrices
     # have a phase freedom and will, therefore, not be equal in general
     a1.deg_shells = None
@@ -49,12 +47,12 @@ with HDFArchive('analyse_block_structure_from_gf.out.h5','r') as ar,\
     assert a1 == a2, 'bs2 not equal'
 
 # check if deg shells are correct
-assert len(SK.deg_shells[0])==1, "wrong number of equivalent groups"
+assert len(SK.deg_shells[0]) == 1, 'wrong number of equivalent groups'
 
 # check if the Green's functions that are found to be equal in the
 # routine are indeed equal
 for d in SK.deg_shells[0]:
-    assert len(d)==2, "wrong number of shells in equivalent group"
+    assert len(d) == 2, 'wrong number of shells in equivalent group'
     # the convention is that for every degenerate shell, the transformation
     # matrix v and the conjugate bool is saved
     # then,
@@ -71,7 +69,7 @@ for d in SK.deg_shells[0]:
         normalized_gfs.append(normalized_gf)
     for i in range(len(normalized_gfs)):
         for j in range(i + 1, len(normalized_gfs)):
-            assert_arrays_are_close(normalized_gfs[i].data, normalized_gfs[j].data, 1.e-5)
+            assert_arrays_are_close(normalized_gfs[i].data, normalized_gfs[j].data, 1.0e-5)
 
 #######################################################################
 # Second test                                                         #
@@ -80,18 +78,21 @@ for d in SK.deg_shells[0]:
 # model                                                               #
 #######################################################################
 
+
 # helper function to get random Hermitian matrix
 def get_random_hermitian(dim):
     herm = np.random.rand(dim, dim) + 1.0j * np.random.rand(dim, dim)
     herm = herm + herm.conjugate().transpose()
     return herm
 
+
 # helper function to get random unitary matrix
 def get_random_transformation(dim):
     herm = get_random_hermitian(dim)
     T = expm(1.0j * herm)
     return T
 
+
 # we will conjugate the Green's function blocks according to the entries
 # of conjugate_values
 # for each of the 5 blocks that will be constructed, there is an entry
@@ -159,23 +160,25 @@ for conjugate in conjugate_values:
             G_noisy['ud'][i : i + 2, i : i + 2] << G_noisy['ud'][i : i + 2, i : i + 2].transpose()
 
     # analyse the block structure
-    G_new = SK.analyse_block_structure_from_gf(G, 1.e-7)
+    G_new = SK.analyse_block_structure_from_gf(G, 1.0e-7)
 
     # transform G_noisy etc. to the new block structure
     G_noisy = SK.block_structure.convert_gf(G_noisy, block_structure1, beta=G_noisy.mesh.beta, space_from='sumk')
     G_pre_transform = SK.block_structure.convert_gf(G_pre_transform, block_structure1, beta=G_noisy.mesh.beta, space_from='sumk')
-    G_noisy_pre_transform = SK.block_structure.convert_gf(G_noisy_pre_transform, block_structure1, beta = G_noisy.mesh.beta, space_from='sumk')
+    G_noisy_pre_transform = SK.block_structure.convert_gf(
+        G_noisy_pre_transform, block_structure1, beta=G_noisy.mesh.beta, space_from='sumk'
+    )
 
-    assert len(SK.deg_shells[0]) == 2, "wrong number of equivalent groups found"
+    assert len(SK.deg_shells[0]) == 2, 'wrong number of equivalent groups found'
-    assert sorted([len(d) for d in SK.deg_shells[0]]) == [2,3], "wrong number of members in the equivalent groups found"
+    assert sorted([len(d) for d in SK.deg_shells[0]]) == [2, 3], 'wrong number of members in the equivalent groups found'
     for d in SK.deg_shells[0]:
         if len(d) == 2:
-            assert 'ud_3' in d, "shell ud_3 missing"
+            assert 'ud_3' in d, 'shell ud_3 missing'
-            assert 'ud_4' in d, "shell ud_4 missing"
+            assert 'ud_4' in d, 'shell ud_4 missing'
         if len(d) == 3:
-            assert 'ud_0' in d, "shell ud_0 missing"
+            assert 'ud_0' in d, 'shell ud_0 missing'
-            assert 'ud_1' in d, "shell ud_1 missing"
+            assert 'ud_1' in d, 'shell ud_1 missing'
-            assert 'ud_2' in d, "shell ud_2 missing"
+            assert 'ud_2' in d, 'shell ud_2 missing'
 
         # the convention is that for every degenerate shell, the transformation
         # matrix v and the conjugate bool is saved
@@ -195,18 +198,18 @@ for conjugate in conjugate_values:
             for j in range(i + 1, len(normalized_gfs)):
                 # here, we use a threshold that is 1 order of magnitude less strict
                 # because of numerics
-                assert_gfs_are_close(normalized_gfs[i], normalized_gfs[j], 1.e-6)
+                assert_gfs_are_close(normalized_gfs[i], normalized_gfs[j], 1.0e-6)
 
     # now we check symm_deg_gf
     # symmetrizing the GF has is has to leave it unchanged
     G_new_symm = G_new[0].copy()
     SK.symm_deg_gf(G_new_symm, 0)
-    assert_block_gfs_are_close(G_new[0], G_new_symm, 1.e-6)
+    assert_block_gfs_are_close(G_new[0], G_new_symm, 1.0e-6)
 
     # symmetrizing the noisy GF, which was carefully constructed,
     # has to give the same result as G_new[0]
     SK.symm_deg_gf(G_noisy, 0)
-    assert_block_gfs_are_close(G_new[0], G_noisy, 1.e-6)
+    assert_block_gfs_are_close(G_new[0], G_noisy, 1.0e-6)
 
     # check backward compatibility of symm_deg_gf
     # first, construct the old format of the deg shells
@@ -217,9 +220,9 @@ for conjugate in conjugate_values:
     # symmetrizing the GF as is has to leave it unchanged
     G_new_symm << G_pre_transform
     SK.symm_deg_gf(G_new_symm, 0)
-    assert_block_gfs_are_close(G_new_symm, G_pre_transform, 1.e-6)
+    assert_block_gfs_are_close(G_new_symm, G_pre_transform, 1.0e-6)
 
     # symmetrizing the noisy GF pre transform, which was carefully constructed,
     # has to give the same result as G_pre_transform
     SK.symm_deg_gf(G_noisy_pre_transform, 0)
-    assert_block_gfs_are_close(G_noisy_pre_transform, G_pre_transform, 1.e-6)
+    assert_block_gfs_are_close(G_noisy_pre_transform, G_pre_transform, 1.0e-6)

test/python/analyse_block_structure_from_gf2.py

@@ -48,7 +48,7 @@ G['ud'] << inverse(inverse(G['ud']) - Hloc)
 
 
 SK = SumkDFT(hdf_file = 'SrIrO3_rot.h5', use_dft_blocks=False)
-G_new = SK.analyse_block_structure_from_gf([G])
+G_new = SK.analyse_block_structure_from_gf([G]*2)
 G_new_symm = G_new[0].copy()
 SK.symm_deg_gf(G_new_symm, 0)
 assert_block_gfs_are_close(G_new[0], G_new_symm)
@@ -93,7 +93,7 @@ known_moments[1,:] = np.eye(10)
 tail, err = fit_tail(G['ud'], known_moments)
 Gt['ud'].set_from_fourier(G['ud'], tail)
 
-G_new = SK.analyse_block_structure_from_gf([Gt])
+G_new = SK.analyse_block_structure_from_gf([Gt] * 2)
 G_new_symm = G_new[0].copy()
 SK.symm_deg_gf(G_new_symm, 0)
 assert_block_gfs_are_close(G_new[0], G_new_symm)

test/python/dc_test.py

@@ -104,15 +104,15 @@ SK_new.set_mu(13.9)
 
 
 icrsh = 0
-dens = SK_compat.density_matrix()
+dens = SK_compat.density_matrix(transform_to_solver_blocks=True)
 
 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_up = np.trace(dens[0]['up_0'].real)
-N_down = np.trace(dens[0]['down'].real)
+N_down = np.trace(dens[0]['down_0'].real)
 N_tot = N_up + N_down
 
 mpi.report(f"{N_up=} ,{N_down=}, {N_tot=}\n")
@@ -141,16 +141,15 @@ 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()
+dens = SK_compat.density_matrix(transform_to_solver_blocks=True)
-
 
 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_up = np.trace(dens[0]['up_0'].real)
-N_down = np.trace(dens[0]['down'].real)
+N_down = np.trace(dens[0]['down_0'].real)
 N_tot = N_up + N_down
 
 mpi.report(f"{N_up=} ,{N_down=}, {N_tot=}\n")

test/python/sumkdft_basic.py

@@ -20,16 +20,20 @@
 #
 ################################################################################
 
-from h5 import *
+from h5 import HDFArchive
 from triqs_dft_tools.sumk_dft_tools import SumkDFTTools
 import triqs.utility.mpi as mpi
 from triqs.utility.comparison_tests import *
 from triqs.utility.h5diff import h5diff
-
+import numpy as np
 SK = SumkDFTTools(hdf_file = 'SrVO3.ref.h5')
 
-dm = SK.density_matrix(method = 'using_gf')
+dm = SK.density_matrix(method = 'using_gf', transform_to_solver_blocks=False, with_Sigma=False)
 dm_pc = SK.partial_charges(with_Sigma=False, with_dc=False)
+dm_pi = SK.density_matrix_using_point_integration()
+
+for key, value in dm[0].items():
+    assert (np.allclose(value, dm_pi[0][key], atol=1e-6, rtol=1e-6))
 
 with HDFArchive('sumkdft_basic.out.h5','w') as ar:
     ar['dm'] = dm