docs: def set_filling(self, populations):
          """Sets the orbital enenergies for on the reference of the free case.
            By setting the desired local populations on every orbital.
            Then generate the necesary operators to respect such configuraion"""
          populations = np.asarray(populations)
  #
  #        self.param['orbital_e'] -= bethe_findfill_zeroT( \
  #                                        self.param['avg_particles'],
  #                                        self.param['orbital_e'],
  #                                        self.param['hopping'])
          efermi = - bethe_find_crystalfield(
                              populations, self.param['hopping'])
          self.param['populations'] = populations
  #        fermion_avg(efermi, self.param['hopping'], 'ocupation')
          self.param['ekin'] = fermion_avg(efermi, self.param['hopping'], 'ekin')
  
          spin_gen_op(self.oper, estimate_gauge(populations))