[cctbxbb] Finding neighbouring atoms using pair_asu_table

[Lukas Reck]

Dear all,

please forgive me if the following is trivial or if I've overlooked
something simple, I'm new to cctbx.

I'm trying to use the cctbx python libraries to read in a large number
of crystal structures and extract some geometric information on the
coordination environment of certain atom types in them. I've been
following the instructions from the 8/2004 CompComm Newsletter No. 4
at <http://www.iucr.org/__data/assets/pdf_file/0003/6384/iucrcompcomm_aug2004.pdf>,
pages 22-28 (is that still up to date?) but I cannot seem to get at
the correct symmetry elements for each neighbouring atom in the
coordination environment of some atom of interest.

As far as I understand, the following should be correct:

assert(isinstance(structure, cctbx.xray.structure)) # yes, it is a structure

asu_mappings = structure.asu_mappings(buffer_thickness = 3) # create
the ASU mappings

pair_asu_table = cctbx.crystal.pair_asu_table(asu_mappings = asu_mappings)
pair_asu_table.add_all_pairs(distance_cutoff = 3) # create the pair ASU table

raw_mappings = asu_mappings.mappings()
raw_pair_table = pair_asu_table.table() # access to the raw structures

for scatterer_index, scatterer in enumerate(structure.scatterers()): #
iterate over all atoms in the structure

    label = scatterer.element_symbol()
    if label == atom_of_interest: # we are only interested in one atom type

        for j_seq, j_syms in raw_pair_table[scatterer_index].items():
# j_seq is the neighbouring scatterer's index, j_sym is a list of
lists of symmetry indices

            for sym_group in j_syms:

                for j_sym in sym_group: # j_sym is an index into
raw_mappings[j_seq]

                    print "j_seq:", j_seq, "j_sym:", j_sym
                    print "i_sym_op:", raw_mappings[j_seq][j_sym].i_sym_op()
                    symmetry = asu_mappings.get_rt_mx(j_seq, j_sym)
                    print "%s at %s" % (scatterer.element_symbol(),
scatterer.site)
                    print "Neighbour: %s at %s" %
(structure.scatterers()[j_seq].element_symbol(),
structure.scatterers()[j_seq].site)
                    print "To be transformed by %s (symmetry #%s)" %
(symmetry, j_sym)
                    print "into %s" %
(symmetry(structure.scatterers()[j_seq].site),)

Here's a sample output for one pair of neighbouring atoms:

j_seq: 10 j_sym: 1
i_sym_op: 10
V at (0.30204, 0.52919, 0.5)
Neighbour: O at (0.36242, 0.49166, 0.5)
To be transformed by y-1/2,-x+1/2,z-1/2 (symmetry #1)
into (-0.00834, 0.13758, 0.0)

As you can see from the coordinates, the neighbouring atom was
correctly identified, but the corresponding symmetry as returned by
get_rt_mx transforms it into a symmetry equivalent that is very far
away.

Clearly I misunderstand the function of one or more of the various
indices here, but I am not quite clear on where my error lies...
Thanks for any pointers!

Best regards

Lukas Reck

--
School of Chemistry
Faculty of Engineering, Mathematics and Science
Trinity College Dublin
Dublin 2

Tel: (+353) 1 896 3452