def to_molden(cartesian_list, buf=None, sort_index=True,
              overwrite=True, float_format='{:.6f}'.format):
    if sort_index:
        cartesian_list = [molecule.sort_index() for molecule in cartesian_list]
    give_header = ("[MOLDEN FORMAT]\n"
                   + "[N_GEO]\n"
                   + str(len(cartesian_list)) + "\n"
                   + '[GEOCONV]\n'
                   + 'energy\n{energy}'
                   + 'max-force\n{max_force}'
                   + 'rms-force\n{rms_force}'
                   + '[GEOMETRIES] (XYZ)\n').format
    values = len(cartesian_list) * '1\n'
    energy = [str(m.metadata.get('energy', 1)) for m in cartesian_list]
    energy = '\n'.join(energy) + '\n'
    header = give_header(energy=energy, max_force=values, rms_force=values)
    coordinates = [x.to_xyz(sort_index=sort_index, float_format=float_format)
                   for x in cartesian_list]
    output = header + '\n'.join(coordinates)
    if buf is not None:
        if overwrite:
            with open(buf, mode='w') as f:
                f.write(output)
        else:
            with open(buf, mode='x') as f:
                f.write(output)
    else:
        return output