code: 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