On Mon, Oct 15, 2012 at 12:51 PM, Dmytro Guzenko
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.)