diff --git a/python/edi2py.py b/python/edi2py.py
index 0e83d826..db3034ef 100644
--- a/python/edi2py.py
+++ b/python/edi2py.py
@@ -11,6 +11,7 @@
 class Link:
     def __init__(self,library):
         self.library = library
+        self.Nineq = None
 
 #function that will add a variable to the dummy class, will be called in variable definition
 def add_global_variable(obj, dynamic_name, target_object, target_attribute):
@@ -130,6 +131,7 @@ def setter(self, new_value):
 import func_io
 global_env.get_sigma = types.MethodType(func_io.get_sigma, global_env)
 global_env.get_gimp = types.MethodType(func_io.get_gimp, global_env)
+global_env.get_dens = types.MethodType(func_io.get_dens, global_env)
 
 #bath_fit
 import func_bath_fit
diff --git a/python/func_io.py b/python/func_io.py
index cd91d2ec..91602317 100644
--- a/python/func_io.py
+++ b/python/func_io.py
@@ -120,4 +120,46 @@ def get_gimp(self,gimp,Nlat=-1,axis="m",typ="n"):
             ed_get_gimp_site_n6(gimp,Dimgimp,Nlat,c_char_p(axis.encode()),c_char_p(typ.encode()))
         else:
             raise ValueError('Shape(array) != 3 in get_bath_component')
-    return gimp
\ No newline at end of file
+    return gimp
+    
+    
+#observables
+
+#density
+def get_dens(self,ilat=None,iorb=None):
+
+    aux_norb = c_int.in_dll(self.library, "Norb").value
+    
+    ed_get_dens_n1_wrap = self.library.ed_get_dens_n1
+    ed_get_dens_n1_wrap.argtypes = [np.ctypeslib.ndpointer(dtype=float,ndim=1, flags='F_CONTIGUOUS')]
+    ed_get_dens_n1_wrap.restype = None
+    
+    ed_get_dens_n2_wrap = self.library.ed_get_dens_n2
+    ed_get_dens_n2_wrap.argtypes = [np.ctypeslib.ndpointer(dtype=float,ndim=2, flags='F_CONTIGUOUS'),
+                                    np.ctypeslib.ndpointer(dtype=np.int64,ndim=1, flags='F_CONTIGUOUS'),
+                                    c_int]
+    ed_get_dens_n2_wrap.restype = None
+
+    if self.Nineq is None:
+        densvec = np.zeros(aux_norb,dtype=float,order="F")
+        ed_get_dens_n1_wrap(densvec)
+        
+        if ilat is not None:
+            raise ValueError("ilat cannot be none for single-impurity DMFT")
+        elif iorb is not None:
+            return densvec[iorb]
+        else:
+            return densvec
+    else:
+        densvec = np.zeros([self.Nineq,aux_norb],dtype=float,order="F")
+        dim_densvector = np.asarray(np.shape(densvec),dtype=np.int64,order="F")
+        ed_get_dens_n2_wrap(densvec,self.Nineq)
+        
+        if ilat is not None and iorb is not None:
+            return densvec[ilat,iorb]
+        elif ilat is None and iorb is not None:
+            return densvec[:,iorb]
+        elif ilat is not None and iorb is None:
+            return densvec[ilat,:]
+        else:
+            return densvec
\ No newline at end of file
diff --git a/python/func_main.py b/python/func_main.py
index 28bc6239..c28bbbe4 100644
--- a/python/func_main.py
+++ b/python/func_main.py
@@ -20,6 +20,7 @@ def init_solver(self,bath):
         init_solver_site(bath,dim_bath)
     else:
         init_solver_ineq(bath,dim_bath)
+        self.Nineq = np.shape(bath)[0] #save number of inequivalent sites: this is useful when we want to know if we are in lattice case or not
 
 # Define the function signature for the Fortran function `solve_site`.
 def solve(self,bath,sflag=True,iflag=True,fmpi=True,mpi_lanc=False):
diff --git a/src/edi2py/edi2py_io.f90 b/src/edi2py/edi2py_io.f90
index 12d7e5f6..6da76df0 100644
--- a/src/edi2py/edi2py_io.f90
+++ b/src/edi2py/edi2py_io.f90
@@ -107,3 +107,20 @@ subroutine ed_get_gimp_lattice_n6_c(self,nlat,d,axis,typ) bind(c, name='ed_get_g
   axis_(1:1)=axis(1)
   call ed_get_gimp(self,nlat,axis_,typ_)
 end subroutine ed_get_gimp_lattice_n6_c
+
+
+!OBSERVABLES
+
+!density
+subroutine ed_get_dens_n1_c(self) bind(c,name="ed_get_dens_n1")
+  real(c_double)     :: self(Norb)
+  call ed_get_dens(self)
+end subroutine ed_get_dens_n1_c
+
+
+subroutine ed_get_dens_n2(self,d,Nlat) bind(c,name="ed_get_dens_n2")
+  integer(c_int64_t) :: d(2)
+  real(c_double)     :: self(d(1),d(2))
+  integer(c_int)     :: Nlat
+  call ed_get_dens(self,Nlat)
+end subroutine ed_get_dens_n2
\ No newline at end of file