The aim of the Xmap is to provide a fast and efficient way of storing some property, in such a way that the crystallographic symmetry and cell repeat are imposed without any effort from the programmer. Therefore, a clipper::Xmap appears to be infinite in every direction, allowing the map to be queried at any position in crystal space. However, only a unique subregion of a single cell is stored. If a value in the map is changed, then every copy of that value throughout crystal space also changes.
There are generally some restrictions on the grid sampling, imposed by the cell dimensions, spacegroup, and FFT requirements. The clipper::Grid_sampling class can be used to select an appropriate grid sampling for any particular problem.
The following image shows a slice of a unit cell which has been sampled on a (6x8) grid. 48 points are stored, with fractional coordinates 0...5/6 on the horizontal axis and 0...7/8 on the vertical axis.
The fraction coordinate of any map grid point may therefore be determined by dividing the zero-based grid coordinate by the grid sampling along each direction.
A naive approach to this problem would be to store only an oblong 'brick' within the unit cell. Unfortunately in the general case this does not work, since some symmetry operators are not aligned along grid axes. Even simple cases, such as a single 2-fold rotation axis, can cause problems, as illustrated in the figure below:
The map section shows a sampling of a P2 unit cell. An attempt has been made to select an asymmetric unit using a rectangular subregion of the section. However it can be seen that the origin is related to itself, and the 2-fold axis also relates the points on the left and right-hand edges of the oblong to other points on the same edge. Therefore 26 points are required to store a complete asymmetric unit section, and the brick containing those points is just over half the unit cell wide and contains 32 points.
The approach used by the Xmap class is therefore to store a compact brick of density, along with a flag array marking which points within that brick are considered to be part of the ASU, and which are considered to be outside the ASU.
Grid coordinates are always integers, the others are floating point values.
Conversion between orthogonal and fractional coordinates is performed using the clipper::Coord_orth::coord_frac() and clipper::Coord_frac::coord_orth() methods, supplying the cell as an argument. e.g. corth.coord_frac(xmap.cell())
Conversion between fractional and grid coordinates is performed using the clipper::Coord_map::coord_frac() and clipper::Coord_frac::coord_map() methods, supplying the grid as an argument. e.g. cfrac.coord_map(xmap.grid_sampling())
The Xmap is constructed by providing a spacegroup, cell, and grid sampling. These are used to determine an appropriate ASU, and allocate an array to store the map data. The data may then be read or written by providing the appropriate grid coordinate. Alternatively, interpolated values and gradients at non-grid locations may be read by providing a fraction coordinate. A variety of interpolation methods are provided for this purpose.
In order to use the class efficiently, some important difficulties must be borne in mind.
A problem arises when we wish to apply some transformation to the values stored in the map. In this case, we must access every unique value in the asymmetric unit once and once only, applying the desired transformation. Only then will the entire map have been transformed correctly.
A second problem arises if we want to access the stored value of the density at some position in crystal space, represented by a grid coordinate. Then it is necessary to search through all the symmetry operators, applying each one in turn to the coordinate to find the operator which brings the coordinate into the stored asymmetric unit. Any cell translations must also be taken into account. The value of the map for that coordinate can then be returned. Clearly this can be time consuming, especially if there are many symmetry operators.
Both these problems are addressed by the use of map reference types. These come in two forms:
The index-like reference behaves like an index: it stores a reference to a map and an index into that map. It is used to loop over all the values in the asymmetric unit of a map, using the Xmap<>::first(), and Map_reference_index::last() and Map_reference_index::next() methods. The coordinate corresponding to the index can be returned at any point.
The coordinate-like reference behaves like a coordinate: it stores a reference to a map and a grid coordinate into that map. However to enhance performance it also stores the index corresponding to that coordinate, and the number of the symmetry operator used to get back into the stored asymmetric unit. Since maps are usually accessed systematically, the next coordinate used will commonly require the same symmetry operator, and so that operator is tried first. An efficient caching mechanism makes incrementing or decrementing the coordinate along the u, v, or w directions particularly fast.
The differences between the index-like and coordinate-like reference types can be summarised as follows:
Use of map reference types is always preferred over requesting a coordinate directly. The accessor methods for the map are designed to encourage efficient usage.
Map reference types may be shared between any maps which have the same symmetry and grid sampling. It is the responsibility of the programmer to ensure this restriction is obeyed.
clipper::CCP4MAPfile file;
file.open_read( "my.map" );
file.import_xmap( xmap );
file.close_read();
The extent and axis order of the input map do not matter, the resulting map will always cover the preferred ASU. However, if there is insufficient information in the input map to obtain an ASU, the remaining portion of the map will be untouched.
To export a map, use:
clipper::CCP4MAPfile file;
file.open_write( "my.map" );
file.export_xmap( xmap );
file.close_write();
The new map is constructed to share a grid and cell with the old map, but with the new spacegroup. We must ensure that every value in the new map is set, so we loop over the new map using a clipper::Xmap_base::Map_reference_index. The we request the corresponding density from the old map, using the coordinate of the Map_reference_index:
clipper::Xmap<float> oldmap, newmap; ... newmap.init( newspacegroup, oldmap.cell(), oldmap.grid_sampling() ); clipper::Xmap_base::Map_reference_index ix; for ( ix = newmap.first(); !ix.last(); ix.next() ) { newmap[ ix ] = oldmap.get_data( ix.coord() ); }
This works, but it would be faster to use a Map_reference_coord to access the second map, so that the symmetry operators do not need to be searched every time. The coordinate of the Map_reference_coord can be set from the Map_reference_index into the first map. i.e.:
clipper::Xmap<float> oldmap, newmap; ... newmap.init( newspacegroup, oldmap.cell(), oldmap.grid_sampling() ); clipper::Xmap_base::Map_reference_index ix; clipper::Xmap_base::Map_reference_coord iy(oldmap); for ( ix = newmap.first(); !ix.last(); ix.next() ) { iy.set_coord( ix.coord() ); newmap[ ix ] = oldmap[ iy ]; }
Three map references can be used, one each for looping along the u, v, and w directions. Suppose the loop must run over the box bounded by Coord_grid g0 and Coord_grid g1. This leads to code of the following form:
clipper::Xmap_base::Map_reference_coord i0, iu, iv, iw; i0 = clipper::Xmap_base::Map_reference_coord( xmap, g0 ); for ( iu = i0; iu.coord().u() <= g1.u(); iu.next_u() ) for ( iv = iu; iv.coord().v() <= g1.v(); iv.next_v() ) for ( iw = iv; iw.coord().w() <= g1.w(); iw.next_w() ) { // ---- access xmap[iw] here ---- }
Alternatively, in most cases it is only necessary to optimise the innermost loop.
int u, v; clipper::Xmap_base::Map_reference_coord ix( xmap ); for ( u = g0.u(); u <= g1.u(); u++ ) for ( v = g0.v(); v <= g1.v(); v++ ) for ( ix.set_coord(Coord_grid(u,v,g0.w())); ix.coord().w() <= g1.w(); ix.next_w() ) { // ---- access xmap[ix] here ---- }