code: def get_plot(self, normalize_rxn_coordinate=True, label_barrier=True):
          
          plt = pretty_plot(12, 8)
          scale = 1 if not normalize_rxn_coordinate else 1 / self.r[-1]
          x = np.arange(0, np.max(self.r), 0.01)
          y = self.spline(x) * 1000
          relative_energies = self.energies - self.energies[0]
          plt.plot(self.r * scale, relative_energies * 1000, 'ro',
                   x * scale, y, 'k-', linewidth=2, markersize=10)
          plt.xlabel("Reaction coordinate")
          plt.ylabel("Energy (meV)")
          plt.ylim((np.min(y) - 10, np.max(y) * 1.02 + 20))
          if label_barrier:
              data = zip(x * scale, y)
              barrier = max(data, key=lambda d: d[1])
              plt.plot([0, barrier[0]], [barrier[1], barrier[1]], 'k--')
              plt.annotate('%.0f meV' % (np.max(y) - np.min(y)),
                           xy=(barrier[0] / 2, barrier[1] * 1.02),
                           xytext=(barrier[0] / 2, barrier[1] * 1.02),
                           horizontalalignment='center')
          plt.tight_layout()
          return plt
docs: Returns the NEB plot. Uses Henkelman's approach of spline fitting
          each section of the reaction path based on tangent force and energies.
  
          Args:
              normalize_rxn_coordinate (bool): Whether to normalize the
                  reaction coordinate to between 0 and 1. Defaults to True.
              label_barrier (bool): Whether to label the maximum barrier.
  
          Returns:
              matplotlib.pyplot object.