The script is intended to be used from PyMol molecular visualization system and allows one to generate &QM_KIND, &LINK and &CELL parts of CP2K input easily. Given a molecular object name and QM part selection name (which can be defined using either approach available in PyMol) it generates output file with corresponding input groups.
To use the script:
1. Save it to some file
2. Start PyMol and issue the following command from it’s shell:
> run /path/to/script.py
3. Select qm part of your system using either mechanism available ni pymol, the issue a command:
> gen_cp2k_qmmm filename, qm_part_selection_name, the_whole_system_object
The output will be saved to file «filename».
N.B.: The file will be overwritten!

The script:

# PyMol script to generate part of the &QMMM section of CP2K input.
# 
# USAGE:
#   gen_cp2k_qmmm filename, qm_part_selection_name, the_whole_system_object
#
#   N.B.: The output is written to file "filename". "Filename" is overwritten!
#
# REFERENCES:
#   CP2K: http://cp2k.berlios.de/
#   PyMol: http://www.pymol.org/
#
# AUTHOR: piton_at_erg.biophys.msu.ru

from pymol import cmd, stored

def gen_cp2k_qmmm (filename, qmsele, mmsele):
  f = open(filename, 'w')
  print >> f, "&QMMM"

# QMMM CELL setup
# n.b. 6A padding
# n.b. cell should be cubic in case of PSOLVER=WAVELET
  ext = cmd.get_extent(qmsele)
  cell = (ext[1][0]-ext[0][0] + 6, ext[1][1]-ext[0][1] + 6, ext[1][2]-ext[0][2] + 6)
  print >> f, "  &CELL\n    ABC %.2f %.2f %.2f\n  &END CELL" % cell

# KINDS
  stored.elems=[]
  cmd.iterate_state(1, qmsele, "stored.elems.append(elem)")
  kinds = set(stored.elems)
  for kind in kinds:
    print >> f, "  &QM_KIND", kind
    stored.ids=[]
    cmd.iterate_state(1, qmsele + ' & elem '+kind, "stored.ids.append(ID)")
    print >> f, "    MM_INDEX", 
    for id in stored.ids: 
      print >> f, id,
    print >> f, "\n  &END QM_KIND"

# LINKS
  stored.qm_ids = []
  cmd.iterate_state(1, qmsele, "stored.qm_ids.append(ID)")
  model = cmd.get_model(mmsele)
  for bond in model.bond:
    [ia, ib] = bond.index
    a = model.atom[ia].id
    b = model.atom[ib].id
    link = 0
    if ((a in stored.qm_ids) and not (b in stored.qm_ids)):
      qmid = a
      mmid = b
      link = 1
    if ((b in stored.qm_ids) and not (a in stored.qm_ids)):
      mmid = a
      qmid = b
      link = 1
    if link: print >> f, "  &LINK\n    QM_INDEX %d\n    MM_INDEX %d\n    QM_KIND H\n  &END LINK" % (qmid, mmid)
 
  print >> f, "&END QMMM"

cmd.extend("gen_cp2k_qmmm", gen_cp2k_qmmm)