Resolve bug in bar_model.delta_f_ creation (alchemistry#397)

* fix bug introduced in PR alchemistry#381: there was a change to creating the delta_f_ matrix,
  which resulted in the columns and indices being tuples that were in the wrong order
  for single lambda computations.
* ensure that columns are in the correct order by explicitly sorting
* add a test for the delta_f_ columns

src/alchemlyb/estimators/bar_.py

@@ -104,41 +104,27 @@ def fit(self, u_nk):
             )
             for i in u_nk.columns
         ]
-<<<<<<< HEAD
 
         # Pull lambda states from indices
         states = list(set(x[1:] if len(x[1:]) > 1 else x[1] for x in u_nk.index))
         for state in states:
-=======
-
-        # Pull lambda states from indices
-        states = list(set( x[1:] for x in u_nk.index))
-        for state in states:
             if len(state) == 1:
                 state = state[0]
->>>>>>> 9eb22a1 (Moving average (#381))
             if state not in self._states_:
                 raise ValueError(
                     f"Indexed lambda state, {state}, is not represented in u_nk columns:"
                     f" {self._states_}"
                 )
-<<<<<<< HEAD
+
         states.sort(key=lambda x: self._states_.index(x))
 
-=======
-
->>>>>>> 9eb22a1 (Moving average (#381))
         # Now get free energy differences and their uncertainties for each step
         deltas = np.array([])
         d_deltas = np.array([])
         for k in range(len(N_k) - 1):
             if N_k[k] == 0 or N_k[k + 1] == 0:
                 continue
-<<<<<<< HEAD
 
-=======
-
->>>>>>> 9eb22a1 (Moving average (#381))
             # get us from lambda step k
             uk = groups.get_group(
                 self._states_[k]
@@ -178,11 +164,7 @@ def fit(self, u_nk):
                 "To compute the free energy with BAR, ensure that values in u_nk exist"
                 f" for the columns:\n{states}."
             )
-<<<<<<< HEAD
 
-=======
-
->>>>>>> 9eb22a1 (Moving average (#381))
         # build matrix of deltas between each state
         adelta = np.zeros((len(deltas) + 1, len(deltas) + 1))
         ad_delta = np.zeros_like(adelta)