code: def selfconsistency(self, u_int, J_coup, mean_field_prev=None):
          
          if mean_field_prev is None:
              mean_field_prev = np.array([self.param['ekin']]*2)
          hlog = [mean_field_prev]
  
          self.oper['Hint'] = self.inter_spin_hamiltonian(u_int, J_coup)
          converging = True
          half_fill = (self.param['populations'] == 0.5).all()
          while converging:
              if half_fill:
                  self.update_H(hlog[-1], self.param['lambda'])
              else:
                  res = root(self.restriction, self.param['lambda'], (hlog[-1]))#, method='lm')
                  if not res.success:
                      res.x = res.x * 0.5 + 0.5*self.param['lambda']
                      self.update_H(self.mean_field()*0.5 + 0.5*hlog[-1], res.x)
                      print('fail', self.param['populations'][3:5])
                  if (self.quasiparticle_weight() < 0.001).all():
                      return hlog
                  self.param['lambda'] = res.x
  
              hlog.append(self.mean_field())
              converging = (abs(hlog[-1] - hlog[-2]) > self.param['tol']).all() \
                      or (abs(self.restriction(self.param['lambda'], hlog[-1])) > self.param['tol']).all()
  
          return hlog
docs: Iterates over the hamiltonian to get the stable selfcosistent one