diff options
| -rw-r--r-- | CMakeLists.txt | 10 | ||||
| -rw-r--r-- | doc/man/cell_tool.1 | 141 | ||||
| -rw-r--r-- | src/cell_tool.c | 306 |
3 files changed, 457 insertions, 0 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 8888be0b..38b1723b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -360,6 +360,15 @@ target_link_libraries(make_pixelmap ${COMMON_LIBRARIES} ${HDF5_C_LIBRARIES}) list(APPEND CRYSTFEL_EXECUTABLES make_pixelmap) # ---------------------------------------------------------------------- +# cell_tool + +set(CELL_TOOL_SOURCES src/cell_tool.c) +add_executable(cell_tool ${CELL_TOOL_SOURCES}) +target_include_directories(cell_tool PRIVATE ${COMMON_INCLUDES}) +target_link_libraries(cell_tool ${COMMON_LIBRARIES}) +list(APPEND CRYSTFEL_EXECUTABLES cell_tool) + +# ---------------------------------------------------------------------- # Install targets set_target_properties(${CRYSTFEL_EXECUTABLES} @@ -391,6 +400,7 @@ install(FILES doc/man/render_hkl.1 doc/man/whirligig.1 doc/man/make_pixelmap.1 + doc/man/cell_tool.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1 ) diff --git a/doc/man/cell_tool.1 b/doc/man/cell_tool.1 new file mode 100644 index 00000000..eb8a4adf --- /dev/null +++ b/doc/man/cell_tool.1 @@ -0,0 +1,141 @@ +.\" +.\" cell_tool man page +.\" +.\" Copyright © 2015-2019 Deutsches Elektronen-Synchrotron DESY, +.\" a research centre of the Helmholtz Association. +.\" +.\" Part of CrystFEL - crystallography with a FEL +.\" + +.TH CELL_TOOL 1 +.SH NAME +cell_tool \- manipulate unit cells +.SH SYNOPSIS +.PP +\fBcell_tool --find-ambi \fImy_structure.cell \fR[\fB-y \fImypointgroup\fR] [\fB--tolerance=\fItols\fR] +.PP +\fBcell_tool --uncenter \fImy_structure.cell \fR[\fB-o \fIoutput.cell\fR] +.PP +\fBcell_tool --rings \fImy_structure.cell \fR[\fB--highres=\fIangstroms\fR] +.PP +\fBcell_tool --compare-cell \fIreference.cell \fImy_structure.cell \fR[\fB--tolerance=\fItols\fR] +.PP +\fBcell_tool --help\fI + +.SH DESCRIPTION +\fBcell_tool\fR performs various manipulations on unit cells, including generating power ring positions, comparing one unit cell to another, calculating a primitive unit cell from a centered one and searching for indexing ambiguities. +.PP +The unit cell can be given as a CrystFEL unit cell file, or alternatively as a PDB file. + +.SH CALCULATING POWDER RING POSITIONS +.PP +\fBcell_tool --rings \fImy_structure.cell \fR[--highres=\fIangstroms\fR] [-y \fIpg\fR] +.PP +This will generate a list of d-spacings and hkl values for the powder rings given by the unit cell file. Note that screw axis and glide plane absences will not be taken into account, so some rings may be absent depending on the space group. +.PP +If you additionally specify the point group using \fB-y\fR (see 'man crystfel' for how to specify point groups), symmetrically equivalent rings will be combined and multiplicities calculated. + +.SH GENERATING A PRIMITIVE UNIT CELL +.PP +\fBcell_tool --uncenter \fImy_structure.cell \fR[-o \fIoutput.cell\fR] +.PP +This will generate a primitive unit cell representing the same lattice as the input. Add the \fB-o\fR option to write the result to a new unit cell file. +.PP +There are an infinite number of primitive unit cell for any lattice. This program generates only one of them. + +.SH COMPARING UNIT CELLS +.PP +\fBcell_tool --compare-cell \fIreference.cell my_structure.cell \fR[\fB--tolerance=\fItols\fR] +.PP +The program will compare the two cells, and report if \fImy_structure.cell\fR can be made to look similar to \fIreference.cell\fR applying any transformation matrix. +.PP +The tolerance \fItols\fR is given as lengthtol,angtol, in percent and degrees respectively, which will be applied to the real-space unit cell axis lengths and angles. +.PP +Note that the comparison is quite limited in the current version. The transformation matrix applied to \fImy_structure.cell\fR is restricted to positive integer (or zero) components, so if \fImy_structure.cell\fR is smaller than the reference, you might need to try swapping the arguments. Even then, it might miss a relationship if one axis is needs to be (say) half the length while another needs to be double. The centering is also ignored. + +.SH FINDING INDEXING AMBIGUITIES +.PP +\fBcell_tool --find-ambi \fImy_structure.cell \fR[-y \fIpg\fR] [\fB--tolerance=\fItols\fR] +.PP +The program will report all transformation matrices which produce a similar unit cell, to within the specified tolerance. The tolerance \fItols\fR is given as lengthtol,angtol, in percent and degrees respectively, which will be applied to the real-space unit cell axis lengths and angles. +.PP +If you additionally give the true symmetry using \fB-y\fR, the program will calculate the ambiguity operators, i.e. the operations which are not symmetry operators of the structure, but which nevertheless leave the lattice looking the same. +.PP +\fBExample 1: Merohedral indexing ambiguity in photosystem I\fR + +The space group of photosystem I crystals as described by PDB code 1JB0 is P63, +so the point group is '6': + +$ cell_tool --find-ambi 1JB0.pdb -y 6 +.nf +[...] +Observed symmetry operations: + 1 : hkl -h-k,k,-l -h-k,h,l -h,-k,l -h,h+k,-l + -k,-h,-l -k,h+k,l k,-h-k,l k,h,-l h,-h-k,-l + h+k,-h,l h+k,-k,-l +Ambiguity operations: + 1 -> 6 : hkl -h-k,k,-l +.fi + +There are 12 reflections which cannot be distinguished between by the lattice alone, but only 6 of those are true symmetry equivalents according to the structure. The transformation describing the indexing ambiguity as follows: "A reflection hkl will be confused with one with indices -h-k,k,-l". Had the point group of the crystals been '622', there would have been no indexing ambiguity (try it!). + +.PP +\fBExample 2: No indexing ambiguity in lysozyme\fR + +The space group of lysozyme crystals as described by PDB code 1VDS is P 43 21 2, so the point group is '422': + +.nf +$ cell_tool --find-ambi 1VDS.pdb -y 422 +[...] +Observed symmetry operations: + 1 : hkl -h,-k,l -h,k,-l -k,-h,-l -k,h,l + k,-h,l k,h,-l h,-k,-l +Ambiguity operations: + 1 -> 422 : hkl +.fi + +All of the apparently equivalent reflections are true symmetry equivalents according to point group 422, so there is no indexing ambiguity. The only operation produced by left coset decomposition is the identity operation. + +.PP +\fBExample 3: "Accidental" ambiguity in myoglobin\fR + +The space group of myoglobin crystals as described by PDB code 3VAU is P2, so the point group is '2'. Note that the "unique axis b" convention has been used for 3VAU (see "man crystfel" for information about specifying point groups): + +.nf +$ cell_tool --find-ambi 3VAU.pdb -y 2_uab +[...] + a b c alpha beta gamma + 3.51 2.84 6.29 nm 90.00 105.50 90.00 deg +[...] + a b c alpha beta gamma + 3.51 2.84 6.33 nm 90.00 106.84 90.00 deg +[...] +Observed symmetry operations: + 1 : hkl -h,-k,h+l -h,k,-l h,-k,-h-l +Ambiguity operations: + 1 -> 2 : hkl -h,-k,h+l +.fi + +The transformations '-h,-k,h+l' and 'h,-k,-h-l', which correspond to indexing "diagonally", produce cells which look very similar to the original cell - a difference of only 0.4A and 1.34 degrees. These two transformations are themselves related by a twofold rotation, which is a true symmetry of this crystal structure. There is therefore only one ambiguity transformation. The transformation is strange because it isn't one of the symmetries displayed by a monoclinic lattice in general. This ambiguity has arisen because of of the particular unit cell parameters for this structure. + +.SH AUTHOR +This page was written by Thomas White. + +.SH REPORTING BUGS +Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>. + +.SH COPYRIGHT AND DISCLAIMER +Copyright © 2015-2019 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association. +.P +cell-tool, and this manual, are part of CrystFEL. +.P +CrystFEL is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +.P +CrystFEL is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +.P +You should have received a copy of the GNU General Public License along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + +.SH SEE ALSO +.BR crystfel (7), +.BR indexamajig (1), +.BR get_hkl (1) diff --git a/src/cell_tool.c b/src/cell_tool.c new file mode 100644 index 00000000..adf9c2a0 --- /dev/null +++ b/src/cell_tool.c @@ -0,0 +1,306 @@ +/* + * cell_tool.c + * + * Unit cell tool + * + * Copyright © 2018 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2012-2018 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdarg.h> +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <unistd.h> +#include <getopt.h> +#include <assert.h> + +#include "cell.h" +#include "cell-utils.h" +#include "reflist-utils.h" +#include "reflist.h" + + +static void show_help(const char *s) +{ + printf("Syntax: %s [options]\n\n", s); + printf( +"Unit cell manipulation tool.\n" +"\n" +" -h, --help Display this help message.\n" +" -p, --pdb=<file> Get unit cell from <file> (PDB or CrystFEL format).\n" +"\n" +" Actions:\n" +" --find-ambi Find indexing ambiguities for the cell.\n" +" --uncenter Calculate a primitive cell.\n" +" --rings Calculate powder ring positions.\n" +" --compare-cell <file> Compare unit cell with cell from <file>.\n" +"\n" +" -y <pointgroup> Real point group of the structure.\n" +" --tolerance=<tol> Set the tolerances for cell reduction.\n" +" Default: 5,5,5,1.5.\n" +); +} + + +static int cells_are_similar(UnitCell *cell1, UnitCell *cell2, float *tols) +{ + double a1, b1, c1, al1, be1, ga1; + double a2, b2, c2, al2, be2, ga2; + UnitCellTransformation *tfn1, *tfn2; + UnitCell *pcell1, *pcell2; + + /* Compare primitive cells, not centered */ + pcell1 = uncenter_cell(cell1, &tfn1); + pcell2 = uncenter_cell(cell2, &tfn2); + + cell_get_parameters(pcell1, &a1, &b1, &c1, &al1, &be1, &ga1); + cell_get_parameters(pcell2, &a2, &b2, &c2, &al2, &be2, &ga2); + + cell_free(pcell1); + cell_free(pcell2); + + if ( !within_tolerance(a1, a2, tols[0]) ) return 0; + if ( !within_tolerance(b1, b2, tols[1]) ) return 0; + if ( !within_tolerance(c1, c2, tols[2]) ) return 0; + if ( !within_tolerance(al1, al2, tols[3]) ) return 0; + if ( !within_tolerance(be1, be2, tols[3]) ) return 0; + if ( !within_tolerance(ga1, ga2, tols[3]) ) return 0; + + return 1; +} + + +static int find_ambi(UnitCell *cell, SymOpList *sym, float *tols) +{ + SymOpList *amb; + SymOpList *ops; + signed int i[9]; + const int maxorder = 3; + + ops = get_pointgroup("1"); + if ( ops == NULL ) return 1; + set_symmetry_name(ops, "Observed"); + + STATUS("Looking for ambiguities up to 3x each lattice length.\n"); + STATUS("This will take about 30 seconds. Please wait...\n"); + + for ( i[0]=-maxorder; i[0]<=+maxorder; i[0]++ ) { + for ( i[1]=-maxorder; i[1]<=+maxorder; i[1]++ ) { + for ( i[2]=-maxorder; i[2]<=+maxorder; i[2]++ ) { + for ( i[3]=-maxorder; i[3]<=+maxorder; i[3]++ ) { + for ( i[4]=-maxorder; i[4]<=+maxorder; i[4]++ ) { + for ( i[5]=-maxorder; i[5]<=+maxorder; i[5]++ ) { + for ( i[6]=-maxorder; i[6]<=+maxorder; i[6]++ ) { + for ( i[7]=-maxorder; i[7]<=+maxorder; i[7]++ ) { + for ( i[8]=-maxorder; i[8]<=+maxorder; i[8]++ ) { + + UnitCellTransformation *tfn; + UnitCell *nc; + IntegerMatrix *m; + int j, k; + int l = 0; + + m = intmat_new(3, 3); + for ( j=0; j<3; j++ ) { + for ( k=0; k<3; k++ ) { + intmat_set(m, j, k, i[l++]); + } + } + + if ( intmat_det(m) != +1 ) continue; + + tfn = tfn_from_intmat(m); + nc = cell_transform(cell, tfn); + + if ( cells_are_similar(cell, nc, tols) ) { + if ( !intmat_is_identity(m) ) add_symop(ops, m); + STATUS("-----------------------------------------------" + "-------------------------------------------\n"); + cell_print(nc); + intmat_print(m); + } else { + intmat_free(m); + } + + tfn_free(tfn); + cell_free(nc); + + } + } + } + } + } + } + } + } + } + + STATUS("Observed symmetry operations:\n"); + describe_symmetry(ops); + + amb = get_ambiguities(ops, sym); + if ( amb == NULL ) { + STATUS("No ambiguities (or error calculating them)\n"); + } else { + STATUS("Ambiguity operations:\n"); + describe_symmetry(amb); + free_symoplist(amb); + } + + free_symoplist(ops); + + return 0; +} + + +enum { + CT_NOTHING, + CT_FINDAMBI, + CT_UNCENTER, + CT_RINGS, + CT_COMPARE, + CT_CHOICES, +}; + + +int main(int argc, char *argv[]) +{ + int c; + char *cell_file = NULL; + UnitCell *cell; + char *toler = NULL; + float tols[4] = {5.0, 5.0, 5.0, 1.5}; /* a,b,c,angles (%,%,%,deg) */ + char *sym_str = NULL; + SymOpList *sym; + int mode = CT_NOTHING; + char *comparecell = NULL; + + /* Long options */ + const struct option longopts[] = { + {"help", 0, NULL, 'h'}, + {"pdb", 1, NULL, 'p'}, + {"tolerance", 1, NULL, 2}, + + /* Modes of operation */ + {"find-ambi", 0, &mode, CT_FINDAMBI}, + {"uncenter", 0, &mode, CT_UNCENTER}, + {"uncentre", 0, &mode, CT_UNCENTER}, + {"rings", 0, &mode, CT_RINGS}, + {"compare-cell", 1, NULL, 3}, + {"highres", 1, NULL, 5}, + + {0, 0, NULL, 0} + }; + + /* Short options */ + while ((c = getopt_long(argc, argv, "hp:y:", + longopts, NULL)) != -1) { + + switch (c) { + + case 'h' : + show_help(argv[0]); + return 0; + + case 'p' : + cell_file = strdup(optarg); + break; + + case 'y' : + sym_str = strdup(optarg); + break; + + case 2 : + toler = strdup(optarg); + break; + + case 3 : + comparecell = strdup(optarg); + mode = CT_COMPARE; + break; + + case 0 : + break; + + default : + return 1; + + } + + } + + if ( cell_file == NULL ) { + ERROR("You must give a filename for the unit cell PDB file.\n"); + return 1; + } + cell = load_cell_from_file(cell_file); + if ( cell == NULL ) { + ERROR("Failed to load cell from '%s'\n", cell_file); + return 1; + } + free(cell_file); + + if ( toler != NULL ) { + int ttt; + ttt = sscanf(toler, "%f,%f,%f,%f", + &tols[0], &tols[1], &tols[2], &tols[3] ); + if ( ttt != 4 ) { + ERROR("Invalid parameters for '--tolerance'\n"); + return 1; + } + free(toler); + } + + STATUS("Input unit cell:\n"); + cell_print(cell); + + if ( validate_cell(cell) ) { + ERROR("Cell is invalid.\n"); + return 1; + } + + if ( sym_str == NULL ) { + ERROR("Please specify the point group of the structure.\n"); + return 1; + } + + sym = get_pointgroup(sym_str); + if ( sym == NULL ) return 1; + free(sym_str); + + if ( mode == CT_NOTHING ) { + ERROR("Please specify mode of operation (see --help)\n"); + return 1; + } + + if ( mode == CT_FINDAMBI ) return find_ambi(cell, sym, tols); + + /* FIXME: Everything else */ + ERROR("Sorry, this mode of operation is not yet implemented.\n"); + return 1; +} |