[cctbxbb] Patterson map data structure

[Nathaniel Echols]

On Mon, Oct 15, 2012 at 12:51 PM, Dmytro Guzenko
<Dmytro.Guzenko at pharm.kuleuven.be> wrote:
> I am playing with the Patterson function as calculated by cctbx and I am
> confused as to how I can extract and manipulate the data.
>
> More specifically, miller.patterson_map(…).real_map() returns a
> one-dimensional array of values. How do I get the corresponding 3D
> coordinates? I suspect it has to do with the gridding and some kind of fixed
> order enumeration, but I haven’t found an example (the peak_search().sites()
> seems to be relevant, but I would like the whole thing).

The map object does indeed have at its core a one-dimensional
flex.double or flex.complex_double object, but with the addition of a
grid accessor, which allows you to index the array using a tuple.  One
could generate a similar object using this python code:

from scitbx.array_family import flex
map_values = flex.random_double(1000)
map_values.reshape(flex.grid(10,10,10))
print map[5,5,5]

Here 5,5,5 is simply the index of the grid point of interest, which
will presumably correspond to some XYZ value, but which one is
somewhat arbitrary depending on the desired grid spacing and the
preferences of the FFT algorithm.  You can get the size of the grid
from the cctbx.miller.fft_map class returned by
array.patterson_map(...) as fft_map.n_real().  Thus:

patt_fft_map = fobs.patterson_map()
patt_map = patt_fft_map.real_map_unpadded()
n_real = patt_fft_map.n_real()
unit_cell = patt_fft_map.unit_cell()

Given a grid coordinate u,v,w you can calculate the fractional
coordinate by dividing u,v,w by n_real, and if necessary, the
Cartesian coordinate by unit_cell.orthogonalize(site_frac=site_frac).
In practice, much of what we do is starting from known 3D coordinates
(fractional or Cartesian) and either pulling out map values using
interpolation, or selecting a set of 1D grid indices around those
coordinates, rather than looking at grid values individually.

I hope this is at least a semi-coherent explanation - if not hopefully
the other developers can correct/clarify.  (I just use this code,
almost none of it is mine!)

-Nat

PS. FYI, the computational overhead for a Python call like
"map[5,5,5]" is relatively large, so if you're looping over 3D
indices, you're probably going to want to use C++ eventually.  (Python
is fine for prototyping, of course.)