from __future__ import division
from cctbx import crystal
import iotbx.pdb
import sys

def run(args):
  assert len(args) > 0
  # read file and construct objects
  pdb_inp = iotbx.pdb.input(file_name=args[0])
  xrs = pdb_inp.xray_structure_simple()
  pdb_h = pdb_inp.construct_hierarchy()
  atoms = pdb_h.atoms()
  cryst_symmetry = xrs.crystal_symmetry()
  # start actual work
  asu_mappings = xrs.asu_mappings(buffer_thickness=5)
  # get all atom pairs within distance_cutoff distance
  pair_generator = crystal.neighbors_fast_pair_generator(
      asu_mappings,
      distance_cutoff=4)
  n_contacts = 0
  sym_operations = []
  for pair in pair_generator:
    # obtain rt_mx_ji - symmetry operator that should be applied to j-th atom
    # to transfer it to i-th atom
    rt_mx_i = asu_mappings.get_rt_mx_i(pair)
    rt_mx_j = asu_mappings.get_rt_mx_j(pair)
    rt_mx_ji = rt_mx_i.inverse().multiply(rt_mx_j)
    # if it is not a unit matrix, that is symmetry related pair of atoms
    if not rt_mx_ji.is_unit_mx():
      if rt_mx_ji not in sym_operations:
        sym_operations.append(rt_mx_ji)
      print pair.i_seq, pair.j_seq, rt_mx_ji, atoms[pair.i_seq].id_str(), 
      print atoms[pair.j_seq].id_str(), "dist=",pair.dist_sq**.5
      n_contacts += 1
  print n_contacts
  for m in sym_operations:
    print m

if (__name__ == "__main__"):
  run(sys.argv[1:])