def main(path_output, s1, s2, dt):
    # Check if the file with couplings exists
    if not os.path.isfile('couplings.npy'):
        # Check all the files stored
        files_im = glob.glob('Ham_*_im')
        # Read the couplings
        couplings = np.stack(
           np.loadtxt(f'Ham_{f}_im') for f in range(len(files_im)))
        # Save the file for fast reading afterwards
        np.save('couplings', couplings)
    else:
        couplings = np.load('couplings.npy')
        ts = np.arange(couplings.shape[0]) * dt
        plt.plot(ts, couplings[:, s1, s2] * r2meV)
        plt.xlabel('Time (fs)')
        plt.ylabel('Energy (meV)')
        plt.show()