Update:

added python interface for get_dens function. It works as a function, so
dens = ed.get_dens() with optional arguments ilat and iorb.

We add a global variable Nineq to the edi2py class, which is set upon
calling init_solver, so we know whether we are in real-space or
single-imp DMFT. This decides which get_dens function to call.

Template for all other get_observable functions.

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

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
+    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

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):

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