def METHOD_NAME(job:ConformersJob, indices=None, temperature=298, unit='kcal/mol', lowest=True):
    molecules, energies, populations = get_main_results(job)
    if isinstance(indices, int):
        N_plot = min(indices, len(energies))
        if lowest:
            indices = list(range(N_plot))
        else:
            indices = np.linspace(0, len(energies)-1, N_plot, dtype=np.int32)
    if indices is None:
        indices = list(range(min(3, len(energies))))
    fig, axes = plt.subplots(1, len(indices), figsize=(12,3))
    if len(indices) == 1:
        axes = [axes]
    for ax, i in zip(axes, indices):
        mol = molecules[i]
        E = energies[i]
        population = populations[i]
        plot_molecule(mol, ax=ax)
        ax.set_title(f"#{i+1}\nΔE = {E:.2f} kcal/mol\nPop.: {population:.3f} (T = {temperature} K)")    