docs: def _initializeFreeEnergies(self, verbose=False, method='zeros'):
          """
          Compute an initial guess at the relative free energies.
  
          OPTIONAL ARGUMENTS
          verbose (boolean) - If True, will print debug information (default: False)
          method (string) - Method for initializing guess at free energies.
          'zeros' - all free energies are initially set to zero
          'mean-reduced-potential' - the mean reduced potential is used
  
          """
  
          if (method == 'zeros'):
              # Use zeros for initial free energies.
              if verbose:
                  print("Initializing free energies to zero.")
              self.f_k[:] = 0.0
          elif (method == 'mean-reduced-potential'):
              # Compute initial guess at free energies from the mean reduced
              # potential from each state
              if verbose:
                  print("Initializing free energies with mean reduced potential for each state.")
              means = np.zeros([self.K], float)
              for k in self.states_with_samples:
                  means[k] = self.u_kn[k, 0:self.N_k[k]].mean()
              if (np.max(np.abs(means)) < 0.000001):
                  print("Warning: All mean reduced potentials are close to zero. If you are using energy differences in the u_kln matrix, then the mean reduced potentials will be zero, and this is expected behavoir.")
              self.f_k = means
          elif (method == 'BAR'):
              # For now, make a simple list of those states with samples.
              initialization_order = np.where(self.N_k > 0)[0]
              # Initialize all f_k to zero.
              self.f_k[:] = 0.0
              # Initialize the rest
              for index in range(0, np.size(initialization_order) - 1):
                  k = initialization_order[index]
                  l = initialization_order[index + 1]
                  # forward work
                  # here, we actually need to distinguish which states are which
                  w_F = (
                      self.u_kn[l,self.x_kindices==k] - self.u_kn[k,self.x_kindices==k])
                      #self.u_kln[k, l, 0:self.N_k[k]] - self.u_kln[k, k, 0:self.N_k[k]])
                      # reverse work
                  w_R = (
                      self.u_kn[k,self.x_kindices==l] - self.u_kn[l,self.x_kindices==l])
                      #self.u_kln[l, k, 0:self.N_k[l]] - self.u_kln[l, l, 0:self.N_k[l]])
  
                  if (len(w_F) > 0 and len(w_R) > 0):
                      # BAR solution doesn't need to be incredibly accurate to
                      # kickstart NR.
                      import pymbar.bar
                      self.f_k[l] = self.f_k[k] + pymbar.bar.BAR(
                          w_F, w_R, relative_tolerance=0.000001, verbose=False, compute_uncertainty=False)
                  else:
                      # no states observed, so we don't need to initialize this free energy anyway, as
                      # the solution is noniterative.
                      self.f_k[l] = 0
  
          else:
              # The specified method is not implemented.
              raise ParameterError('Method ' + method + ' unrecognized.')
  
          # Shift all free energies such that f_0 = 0.
          self.f_k[:] = self.f_k[:] - self.f_k[0]
  
          return