diff options
author | Thomas White <taw@physics.org> | 2014-09-19 16:07:24 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2014-09-19 16:23:39 +0200 |
commit | 6a6cb3b4d7f15c234a79ff8421a0ae5c1a1dcb2a (patch) | |
tree | 00f6e0da9a8d086af18b0b1f34433bc115c9f206 | |
parent | 2c959daa7a46b99a10dd5a1998b62ccb8def97de (diff) |
Introduce CrystFEL unit cell files
-rw-r--r-- | Makefile.am | 3 | ||||
-rw-r--r-- | doc/examples/cell-example.cell | 12 | ||||
-rw-r--r-- | doc/man/check_hkl.1 | 3 | ||||
-rw-r--r-- | doc/man/compare_hkl.1 | 3 | ||||
-rw-r--r-- | doc/man/get_hkl.1 | 3 | ||||
-rw-r--r-- | doc/man/indexamajig.1 | 13 | ||||
-rw-r--r-- | doc/man/partial_sim.1 | 9 | ||||
-rw-r--r-- | doc/man/pattern_sim.1 | 5 | ||||
-rw-r--r-- | doc/man/render_hkl.1 | 7 | ||||
-rw-r--r-- | libcrystfel/src/cell-utils.c | 262 | ||||
-rw-r--r-- | libcrystfel/src/cell-utils.h | 9 | ||||
-rw-r--r-- | libcrystfel/src/cell.c | 83 | ||||
-rw-r--r-- | libcrystfel/src/cell.h | 12 | ||||
-rw-r--r-- | libcrystfel/src/dirax.c | 6 | ||||
-rw-r--r-- | libcrystfel/src/mosflm.c | 6 | ||||
-rw-r--r-- | libcrystfel/src/xds.c | 9 | ||||
-rw-r--r-- | src/check_hkl.c | 15 | ||||
-rw-r--r-- | src/compare_hkl.c | 14 | ||||
-rw-r--r-- | src/get_hkl.c | 24 | ||||
-rw-r--r-- | src/indexamajig.c | 16 | ||||
-rw-r--r-- | src/partial_sim.c | 2 | ||||
-rw-r--r-- | src/pattern_sim.c | 4 | ||||
-rw-r--r-- | src/render_hkl.c | 20 |
23 files changed, 420 insertions, 120 deletions
diff --git a/Makefile.am b/Makefile.am index dbb15aa7..f6c13f25 100644 --- a/Makefile.am +++ b/Makefile.am @@ -162,7 +162,8 @@ crystfeldoc_DATA = doc/twin-calculator.pdf doc/examples/lcls-dec.geom \ doc/examples/lcls-xpp-estimate.beam \ doc/examples/simple.geom \ doc/examples/cspad-feb2011.geom \ - doc/examples/lcls-cxi-9keV.beam + doc/examples/lcls-cxi-9keV.beam \ + doc/examples/cell-example.cell EXTRA_DIST += $(crystfeldoc_DATA) doc/twin-calculator.odt \ doc/reference/libcrystfel/xml/coding-standards.xml diff --git a/doc/examples/cell-example.cell b/doc/examples/cell-example.cell new file mode 100644 index 00000000..ac236214 --- /dev/null +++ b/doc/examples/cell-example.cell @@ -0,0 +1,12 @@ +CrystFEL unit cell file version 1.0 + +lattice_type = cubic +centering = I + +a = 66.2 A +b = 66.2 A +c = 66.2 A + +al = 90.0 deg +be = 90.0 deg +ga = 90.0 deg diff --git a/doc/man/check_hkl.1 b/doc/man/check_hkl.1 index 7e21dd10..4853a482 100644 --- a/doc/man/check_hkl.1 +++ b/doc/man/check_hkl.1 @@ -20,10 +20,11 @@ check_hkl calculates figures of merit for reflection data, such as completeness .SH OPTIONS .PD 0 +.IP "\fB-p\fR \fIunitcell.cell\fR" .IP "\fB-p\fR \fIunitcell.pdb\fR" .IP \fB--pdb=\fR\fIunitcell.pdb\fR .PD -Specify the name of the PDB file containing at least a CRYST1 line describing the unit cell. +Specify the name of the file containing unit cell information, in PDB or CrystFEL format. .PD 0 .IP "\fB-y\fR \fIpointgroup\fR" diff --git a/doc/man/compare_hkl.1 b/doc/man/compare_hkl.1 index b2374b73..0ae6a890 100644 --- a/doc/man/compare_hkl.1 +++ b/doc/man/compare_hkl.1 @@ -26,10 +26,11 @@ compare_hkl compares two sets of reflection data and calculates figures of merit Specify the symmetry of the reflections. The symmetry must be the same for both lists of reflections. Default: 1 (no symmetry). .PD 0 +.IP "\fB-p\fR \fIunitcell.cell\fR" .IP "\fB-p\fR \fIunitcell.pdb\fR" .IP \fB--pdb=\fR\fIunitcell.pdb\fR .PD -Specify the name of the PDB file containing at least a CRYST1 line describing the unit cell. +Specify the name of the file containing unit cell information, in PDB or CrystFEL format. .PD 0 .IP \fB--fom=\fR\fIFoM\fR diff --git a/doc/man/get_hkl.1 b/doc/man/get_hkl.1 index 92b8ce37..94c73579 100644 --- a/doc/man/get_hkl.1 +++ b/doc/man/get_hkl.1 @@ -77,10 +77,11 @@ In the first form, reflections with d (=lamba/2*sin(theta)) < \fIn\fR will be re In the second form, anisotropic truncation will be performed with separate resolution limits \fIn1\fR, \fIn2\fR and \fIn3\fR along a*, b* and c* respectively. You must also specify \fB-p\fR or \fB--pdb\fR. .PD 0 +.IP "\fB-p\fR \fIunitcell.cell\fR" .IP "\fB-p\fR \fIunitcell.pdb\fR" .IP \fB--pdb=\fR\fIunitcell.pdb\fR .PD -Specify the name of the PDB file containing at least a CRYST1 line describing the unit cell. +Specify the name of the file containing unit cell information, in PDB or CrystFEL format. .SH AUTHOR This page was written by Thomas White. diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1 index 4b90d762..71f836c6 100644 --- a/doc/man/indexamajig.1 +++ b/doc/man/indexamajig.1 @@ -22,7 +22,7 @@ indexamajig \- bulk indexing and data reduction program \fBindexamajig\fR takes a list of diffraction snapshots from crystals in random orientations and attempts to find peaks, index and integrate each one. The input is a list of diffraction image files in HDF5 format and some auxiliary files and parameters. The output is a long text file ('stream') containing the results from each image in turn. -For minimal basic use, you need to provide the list of diffraction patterns, the method which will be used to index, a file describing the geometry of the detector, and a PDB file which contains the unit cell which will be used for the indexing. Here is what the minimal use might look like on the command line: +For minimal basic use, you need to provide the list of diffraction patterns, the method which will be used to index, a file describing the geometry of the detector, and a file which contains the unit cell which will be used for the indexing. Here is what the minimal use might look like on the command line: .IP \fBindexamajig\fR .PD @@ -57,7 +57,7 @@ To use this option, 'dirax' must be in your shell's search path. If you see the .IP \fBmosflm\fR .PD -As \fBdirax\fR, but invoke MOSFLM instead. If you provide a PDB file (with \fB-p\fR), the lattice type and centering information will be passed to MOSFLM, which will then return solutions which match. Note that the lattice parameter information will \fBnot\fR be given to MOSFLM, because it has no way to make use of it. +As \fBdirax\fR, but invoke MOSFLM instead. If you provide a unit cell (with \fB-p\fR), the lattice type and centering information will be passed to MOSFLM, which will then return solutions which match. Note that the lattice parameter information will \fBnot\fR be given to MOSFLM, because it has no way to make use of it. .sp To use this option, 'ipmosflm' must be in your shell's search path. If you see the MOSFLM version and copyright information when you run \fBipmosflm\fR on the command line, things are set up correctly. @@ -96,7 +96,7 @@ Do not check that the cell accounts for any of the peaks as described in \fBdira .IP \fB-nolatt\fR .PD -Do not use the lattice type information from the PDB file to help guide the indexing. Use with \fBmosflm\fR, which is the only indexing method which can optionally take advantage of this information. This is the default behaviour for \fBdirax\fR. This option makes no sense for \fBreax\fR, which is intrinsically based on using known lattice information. +Do not use the lattice type information to help guide the indexing. Use with \fBmosflm\fR, which is the only indexing method which can optionally take advantage of this information. This is the default behaviour for \fBdirax\fR. This option makes no sense for \fBreax\fR, which is intrinsically based on using known lattice information. .IP \fB-latt\fR .PD @@ -198,10 +198,11 @@ Read the detector geometry description from \fIfilename\fR. See \fBman crystfel Read the beam description from \fIfilename\fR. See \fBman crystfel_geometry\fR for more information. .PD 0 -.IP "\fB-p\fR \fIfilename\fR" -.IP \fB--pdb=\fR\fIfilename\fR +.IP "\fB-p\fR \fIunitcell.cell\fR" +.IP "\fB-p\fR \fIunitcell.pdb\fR" +.IP \fB--pdb=\fR\fIunitcell.pdb\fR .PD -Read the unit cell for comparison from the CRYST1 line of the PDB file called \fIfilename\fR. +Specify the name of the file containing unit cell information, in PDB or CrystFEL format. .PD 0 .IP "\fB-e\fR \fIpath\fR" diff --git a/doc/man/partial_sim.1 b/doc/man/partial_sim.1 index 2af472e9..0c30ba73 100644 --- a/doc/man/partial_sim.1 +++ b/doc/man/partial_sim.1 @@ -25,7 +25,7 @@ partial_sim \- calculate partial reflections partial_sim calculates the intensities of idealised partial reflections from crystals in random orientations, which is useful for testing the convergence of Monte Carlo integration or scaling/post-refinement techniques. .P -You need to provide a CrystFEL geometry file (with \fB--geometry=\fR\fImy.geom\fR or \fB-g\fR \fImy.geom\fR), a beam description file (with \fB--beam=\fR\fImy.beam\fR or \fB-b\fR \fImy.beam\fR), a PDB file containing at least a CRYST1 line specifying the unit cell to use for the simulation (with \fB--pdb=\fR\fImy.pdb\fR or \fB-p\fR \fImy.pdb\fR), and an output filename with \fB--output=\fR\fImy.stream\fR or \fB-o\fR \fImy.stream\fR. +You need to provide a CrystFEL geometry file (with \fB--geometry=\fR\fImy.geom\fR or \fB-g\fR \fImy.geom\fR), a beam description file (with \fB--beam=\fR\fImy.beam\fR or \fB-b\fR \fImy.beam\fR), a file containing the unit cell to use for the simulation (with \fB--pdb=\fR\fImy.pdb\fR or \fB-p\fR \fImy.pdb\fR), and an output filename with \fB--output=\fR\fImy.stream\fR or \fB-o\fR \fImy.stream\fR. For each randomly generated orientation, partial_sim calculates which reflections would appear on the detector with the specified beam parameters. It calculates the partiality for each reflection and multiplies it by the fully integrated intensity to produce a partial intensity. The fully integrated intensities can be taken from a file you provide (see below), otherwise they will be randomly generated (by taking the absolute value of a Gaussian random number, mean zero and standard deviation 1000). All the partial intensities for the orientation are multiplied by an overall scaling factor, which is randomly generated with a Gaussian distribution with mean 1 and standard deviation 0.3. The partial intensities are written to the output stream, and the process repeated for as many different orientations as you ask for (see below, default: 2). @@ -48,6 +48,13 @@ Take the fully integrated reflection intensities from \fIfile.hkl\fR, instead of Specify the number of different orientations to simulate. Default: 2. .PD 0 +.IP "\fB-p\fR \fIunitcell.cell\fR" +.IP "\fB-p\fR \fIunitcell.pdb\fR" +.IP \fB--pdb=\fR\fIunitcell.pdb\fR +.PD +Specify the name of the file containing unit cell information, in PDB or CrystFEL format. + +.PD 0 .B .IP "-r \fIrandom.hkl\fR" .B diff --git a/doc/man/pattern_sim.1 b/doc/man/pattern_sim.1 index e71413db..1c3a3b41 100644 --- a/doc/man/pattern_sim.1 +++ b/doc/man/pattern_sim.1 @@ -25,17 +25,18 @@ pattern_sim simulates diffraction patterns from small crystals probed with femto pattern_sim -g mydetector.geom -b my.beam -p my.pdb -r -i myintensities.hkl -The unit cell geometry will be taken from the CRYST1 line in the PDB file you provide, and the intensities of the reflections will be interpolated from the reflection list file you provide. The reflection list format is the same as that output by process_hkl and handled by get_hkl. You also need beam and geometry description files (-b and -g respectively). See `man crystfel_geometry' for details of how to create a geometry file. Examples of both files can be found in the installation directory, which is normally /usr/local/share/doc/crystfel. +The unit cell geometry will be taken from the unit cell file you provide, and the intensities of the reflections will be interpolated from the reflection list file you provide. The reflection list format is the same as that output by process_hkl and handled by get_hkl. You also need beam and geometry description files (-b and -g respectively). See `man crystfel_geometry' for details of how to create a geometry file. Examples of both files can be found in the installation directory, which is normally /usr/local/share/doc/crystfel. The result will be written to an HDF5 file in the current directory with the name `sim.h5'. .SH OPTIONS .PD 0 +.IP "\fB-p\fR \fIunitcell.cell\fR" .IP "\fB-p\fR \fIunitcell.pdb\fR" .IP \fB--pdb=\fR\fIunitcell.pdb\fR .PD -Specify the name of the PDB file containing at least a CRYST1 line describing the unit cell. +Specify the name of the file containing unit cell information, in PDB or CrystFEL format. .PD 0 .IP \fB--gpu\fR diff --git a/doc/man/render_hkl.1 b/doc/man/render_hkl.1 index ff55e9f0..bd8b40c8 100644 --- a/doc/man/render_hkl.1 +++ b/doc/man/render_hkl.1 @@ -46,10 +46,11 @@ hence represents the Laue zone number. Write the output (in PDF format) to \fIfilename\fR. Default: \fB--output=za.pdf\fR. .PD 0 -.IP "\fB-p\fR \fIfilename\fR" -.IP \fB--pdb=\fR\fIfilenamefR +.IP "\fB-p\fR \fIunitcell.cell\fR" +.IP "\fB-p\fR \fIunitcell.pdb\fR" +.IP \fB--pdb=\fR\fIunitcell.pdb\fR .PD -Get the unit cell parameters from the CRYST1 line contained in \fIfilename\fR. +Specify the name of the file containing unit cell information, in PDB or CrystFEL format. .PD 0 .IP \fB--boost=\fR\fIn\fR diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 32e3bbe7..1bc937f9 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -3,13 +3,13 @@ * * Unit Cell utility functions * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2012 Thomas White <taw@physics.org> - * 2012 Lorenzo Galli + * 2009-2012,2014 Thomas White <taw@physics.org> + * 2012 Lorenzo Galli * * This file is part of CrystFEL. * @@ -128,6 +128,27 @@ LatticeType lattice_from_str(const char *s) } +static int check_centering(char cen) +{ + switch ( cen ) { + + case 'P' : + case 'A' : + case 'B' : + case 'C' : + case 'I' : + case 'F' : + case 'R' : + case 'H' : + return 0; + + default: + return 1; + + } +} + + int right_handed(UnitCell *cell) { double asx, asy, asz; @@ -201,39 +222,46 @@ void cell_print(UnitCell *cell) STATUS(", unique axis %c", cell_get_unique_axis(cell)); } - if ( right_handed(cell) ) { - STATUS(", right handed"); - } else { - STATUS(", left handed"); + if ( cell_has_parameters(cell) ) { + if ( right_handed(cell) ) { + STATUS(", right handed"); + } else { + STATUS(", left handed"); + } } STATUS(", point group '%s'.\n", cell_get_pointgroup(cell)); - cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); + if ( cell_has_parameters(cell) ) { + cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); + + STATUS("a b c alpha beta gamma\n"); + STATUS("%6.2f %6.2f %6.2f A %6.2f %6.2f %6.2f deg\n", + a*1e10, b*1e10, c*1e10, + rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); - STATUS("a b c alpha beta gamma\n"); - STATUS("%6.2f %6.2f %6.2f A %6.2f %6.2f %6.2f deg\n", - a*1e10, b*1e10, c*1e10, - rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + STATUS("a = %10.3e %10.3e %10.3e m\n", ax, ay, az); + STATUS("b = %10.3e %10.3e %10.3e m\n", bx, by, bz); + STATUS("c = %10.3e %10.3e %10.3e m\n", cx, cy, cz); - STATUS("a = %10.3e %10.3e %10.3e m\n", ax, ay, az); - STATUS("b = %10.3e %10.3e %10.3e m\n", bx, by, bz); - STATUS("c = %10.3e %10.3e %10.3e m\n", cx, cy, cz); + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); + STATUS("a* = %10.3e %10.3e %10.3e m^-1 (modulus %10.3e m^-1)\n", + asx, asy, asz, modulus(asx, asy, asz)); + STATUS("b* = %10.3e %10.3e %10.3e m^-1 (modulus %10.3e m^-1)\n", + bsx, bsy, bsz, modulus(bsx, bsy, bsz)); + STATUS("c* = %10.3e %10.3e %10.3e m^-1 (modulus %10.3e m^-1)\n", + csx, csy, csz, modulus(csx, csy, csz)); - STATUS("astar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", - asx, asy, asz, modulus(asx, asy, asz)); - STATUS("bstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", - bsx, bsy, bsz, modulus(bsx, bsy, bsz)); - STATUS("cstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", - csx, csy, csz, modulus(csx, csy, csz)); + STATUS("Cell representation is %s.\n", cell_rep(cell)); - STATUS("Cell representation is %s.\n", cell_rep(cell)); + } else { + STATUS("Unit cell parameters are not specified.\n"); + } } @@ -1044,6 +1072,14 @@ static void determine_lattice(UnitCell *cell, } +/** + * load_cell_from_pdb: + * + * Loads a unit cell from a PDB file. + * + * Returns: a newly allocated %UnitCell. + * + */ UnitCell *load_cell_from_pdb(const char *filename) { FILE *fh; @@ -1116,6 +1152,178 @@ UnitCell *load_cell_from_pdb(const char *filename) } +static int get_length_m(char **bits, int nbits, double *pl) +{ + char *rval; + + if ( nbits < 4 ) { + ERROR("No units specified for '%s'\n", bits[0]); + return 1; + } + + *pl = strtod(bits[2], &rval); + if ( *rval != '\0' ) { + ERROR("Invalid value '%s'.\n", bits[2]); + return 1; + } + + if ( strcmp(bits[3], "nm") == 0 ) { + *pl *= 1e-9; + } else if ( strcmp(bits[3], "A") == 0 ) { + *pl *= 1e-10; + } else { + ERROR("Unrecognised length units '%s'\n", bits[3]); + return 1; + } + + return 0; +} + + +static int get_angle_rad(char **bits, int nbits, double *pl) +{ + char *rval; + + if ( nbits < 4 ) { + ERROR("No units specified for '%s'\n", bits[0]); + return 1; + } + + *pl = strtod(bits[2], &rval); + if ( *rval != '\0' ) { + ERROR("Invalid value '%s'.\n", bits[2]); + return 1; + } + + if ( strcmp(bits[3], "rad") == 0 ) { + /* Do nothing, already in rad */ + } else if ( strcmp(bits[3], "deg") == 0 ) { + *pl = deg2rad(*pl); + } else { + ERROR("Unrecognised angle units '%s'\n", bits[3]); + return 1; + } + + return 0; +} + + +/** + * load_cell_from_file: + * + * Loads a unit cell from a file of any type (PDB or CrystFEL format) + * + * Returns: a newly allocated %UnitCell. + * + */ +UnitCell *load_cell_from_file(const char *filename) +{ + FILE *fh; + char *rval; + char line[1024]; + UnitCell *cell; + int have_a = 0, have_b = 0, have_c = 0; + int have_al = 0, have_be = 0, have_ga = 0; + double a, b, c, al, be, ga; + + fh = fopen(filename, "r"); + if ( fh == NULL ) { + ERROR("Couldn't open '%s'\n", filename); + return NULL; + } + + rval = fgets(line, 1023, fh); + chomp(line); + + if ( strcmp(line, "CrystFEL unit cell file version 1.0") != 0 ) { + fclose(fh); + return load_cell_from_pdb(filename); + } + + cell = cell_new(); + + do { + + char line[1024]; + int n1; + int i; + char **bits; + + rval = fgets(line, 1023, fh); + chomp(line); + + n1 = assplode(line, " \t", &bits, ASSPLODE_NONE); + if ( n1 < 3 ) { + for ( i=0; i<n1; i++ ) free(bits[i]); + free(bits); + continue; + } + + if ( bits[0][0] == ';' ) { + for ( i=0; i<n1; i++ ) free(bits[i]); + free(bits); + continue; + } + + if ( strcmp(bits[0], "lattice_type") == 0 ) { + LatticeType lt = lattice_from_str(bits[2]); + cell_set_lattice_type(cell, lt); + + } else if ( strcmp(bits[0], "centering") == 0 ) { + char cen = bits[2][0]; + if ( !check_centering(cen) ) { + cell_set_centering(cell, cen); + } else { + ERROR("Unrecognised centering '%c'\n", cen); + } + + } else if ( strcmp(bits[0], "a") == 0 ) { + if ( !get_length_m(bits, n1, &a) ) { + have_a = 1; + } + + } else if ( strcmp(bits[0], "b") == 0 ) { + if ( !get_length_m(bits, n1, &b) ) { + have_b = 1; + } + + } else if ( strcmp(bits[0], "c") == 0 ) { + if ( !get_length_m(bits, n1, &c) ) { + have_c = 1; + } + + } else if ( strcmp(bits[0], "al") == 0 ) { + if ( !get_angle_rad(bits, n1, &al) ) { + have_al = 1; + } + + } else if ( strcmp(bits[0], "be") == 0 ) { + if ( !get_angle_rad(bits, n1, &be) ) { + have_be = 1; + } + + } else if ( strcmp(bits[0], "ga") == 0 ) { + if ( !get_angle_rad(bits, n1, &ga) ) { + have_ga = 1; + } + + } else { + ERROR("Unrecognised field '%s'\n", bits[0]); + } + + for ( i=0; i<n1; i++ ) free(bits[i]); + free(bits); + + } while ( rval != NULL ); + + if ( have_a && have_b && have_c && have_al && have_be && have_ga ) { + cell_set_parameters(cell, a, b, c, al, be, ga); + } + + return cell; +} + + /* Force the linker to bring in CBLAS to make GSL happy */ void cell_fudge_gslcblas() { diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index 00071c59..0b900096 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -3,13 +3,13 @@ * * Unit Cell utility functions * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2012 Thomas White <taw@physics.org> - * 2012 Lorenzo Galli + * 2009-2013,2014 Thomas White <taw@physics.org> + * 2012 Lorenzo Galli * * This file is part of CrystFEL. * @@ -56,6 +56,7 @@ extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose, extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *tempcell); extern UnitCell *load_cell_from_pdb(const char *filename); +extern UnitCell *load_cell_from_file(const char *filename); extern int cell_is_sensible(UnitCell *cell); diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 44c767f0..490a80ac 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -3,15 +3,15 @@ * * A class representing a unit cell * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2012 Thomas White <taw@physics.org> - * 2010 Richard Kirian - * 2012 Lorenzo Galli + * 2009-2012,2014 Thomas White <taw@physics.org> + * 2010 Richard Kirian + * 2012 Lorenzo Galli * * This file is part of CrystFEL. * @@ -70,6 +70,11 @@ struct _unitcell { CellRepresentation rep; + int have_parameters; + int have_a; + int have_b; + int have_c; + /* Crystallographic representation */ double a; /* m */ double b; /* m */ @@ -116,9 +121,9 @@ UnitCell *cell_new() cell->a = 1.0; cell->b = 1.0; cell->c = 1.0; - cell->alpha = M_PI_2; - cell->beta = M_PI_2; - cell->gamma = M_PI_2; + cell->alpha = 0.0; + cell->beta = 0.0; + cell->gamma = 0.0; cell->rep = CELL_REP_CRYST; @@ -126,6 +131,10 @@ UnitCell *cell_new() cell->lattice_type = L_TRICLINIC; cell->centering = 'P'; cell->unique_axis = '?'; + cell->have_parameters = 0; + cell->have_a = 0; + cell->have_b = 0; + cell->have_c = 0; return cell; } @@ -146,6 +155,20 @@ void cell_free(UnitCell *cell) } +/** + * cell_has_parameters: + * @cell: A %UnitCell + * + * Returns: True if @cell has its parameters specified. + * + */ +int cell_has_parameters(UnitCell *cell) +{ + if ( cell == NULL ) return 0; + return cell->have_parameters; +} + + void cell_set_parameters(UnitCell *cell, double a, double b, double c, double alpha, double beta, double gamma) { @@ -159,6 +182,7 @@ void cell_set_parameters(UnitCell *cell, double a, double b, double c, cell->gamma = gamma; cell->rep = CELL_REP_CRYST; + cell->have_parameters = 1; } @@ -174,6 +198,7 @@ void cell_set_cartesian(UnitCell *cell, cell->cx = cx; cell->cy = cy; cell->cz = cz; cell->rep = CELL_REP_CART; + cell->have_parameters = 1; } @@ -182,6 +207,10 @@ void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az) if ( cell == NULL ) return; cell->ax = ax; cell->ay = ay; cell->az = az; cell->rep = CELL_REP_CART; + cell->have_a = 1; + if ( cell->have_a && cell->have_b && cell->have_c ) { + cell->have_parameters = 1; + } } @@ -190,6 +219,10 @@ void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz) if ( cell == NULL ) return; cell->bx = bx; cell->by = by; cell->bz = bz; cell->rep = CELL_REP_CART; + cell->have_b = 1; + if ( cell->have_a && cell->have_b && cell->have_c ) { + cell->have_parameters = 1; + } } @@ -198,6 +231,10 @@ void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz) if ( cell == NULL ) return; cell->cx = cx; cell->cy = cy; cell->cz = cz; cell->rep = CELL_REP_CART; + cell->have_c = 1; + if ( cell->have_a && cell->have_b && cell->have_c ) { + cell->have_parameters = 1; + } } @@ -228,6 +265,7 @@ UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs, cell->cxs = cs.u; cell->cys = cs.v; cell->czs = cs.w; cell->rep = CELL_REP_RECIP; + cell->have_parameters = 1; return cell; } @@ -245,6 +283,7 @@ UnitCell *cell_new_from_direct_axes(struct rvec a, struct rvec b, struct rvec c) cell->cx = c.u; cell->cy = c.v; cell->cz = c.w; cell->rep = CELL_REP_CART; + cell->have_parameters = 1; return cell; } @@ -280,6 +319,7 @@ void cell_set_reciprocal(UnitCell *cell, cell->cxs = csx; cell->cys = csy; cell->czs = csz; cell->rep = CELL_REP_RECIP; + cell->have_parameters = 1; } @@ -317,6 +357,11 @@ static int cell_crystallographic_to_cartesian(UnitCell *cell, { double tmp, V, cosalphastar, cstar; + if ( !cell->have_parameters ) { + ERROR("Unit cell has unspecified parameters.\n"); + return 1; + } + /* Firstly: Get a in terms of x, y and z * +a (cryst) is defined to lie along +x (cart) */ *ax = cell->a; @@ -434,6 +479,11 @@ int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c, if ( cell == NULL ) return 1; + if ( !cell->have_parameters ) { + ERROR("Unit cell has unspecified parameters.\n"); + return 1; + } + switch ( cell->rep ) { case CELL_REP_CRYST: @@ -490,6 +540,11 @@ int cell_get_cartesian(UnitCell *cell, { if ( cell == NULL ) return 1; + if ( !cell->have_parameters ) { + ERROR("Unit cell has unspecified parameters.\n"); + return 1; + } + switch ( cell->rep ) { case CELL_REP_CRYST: @@ -526,8 +581,14 @@ int cell_get_reciprocal(UnitCell *cell, { int r; double ax, ay, az, bx, by, bz, cx, cy, cz; + if ( cell == NULL ) return 1; + if ( !cell->have_parameters ) { + ERROR("Unit cell has unspecified parameters.\n"); + return 1; + } + switch ( cell->rep ) { case CELL_REP_CRYST: @@ -697,9 +758,9 @@ UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t) out = cell_new_from_cell(cell); if ( out == NULL ) return NULL; - cell_get_cartesian(out, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz); + if ( cell_get_cartesian(out, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz) ) return NULL; m = gsl_matrix_alloc(3,3); a = gsl_matrix_calloc(3,3); diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index 0a512a43..eecf6007 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -3,15 +3,15 @@ * * A class representing a unit cell * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2012 Thomas White <taw@physics.org> - * 2010,2012 Richard Kirian - * 2012 Lorenzo Galli + * 2009-2012,2014 Thomas White <taw@physics.org> + * 2010,2012 Richard Kirian + * 2012 Lorenzo Galli * * This file is part of CrystFEL. * @@ -117,6 +117,8 @@ extern UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs, extern UnitCell *cell_new_from_direct_axes(struct rvec as, struct rvec bs, struct rvec cs); +extern int cell_has_parameters(UnitCell *cell); + extern void cell_set_cartesian(UnitCell *cell, double ax, double ay, double az, double bx, double by, double bz, diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c index dbd43267..e61a4bfb 100644 --- a/libcrystfel/src/dirax.c +++ b/libcrystfel/src/dirax.c @@ -638,9 +638,9 @@ IndexingPrivate *dirax_prepare(IndexingMethod *indm, UnitCell *cell, if ( *indm & INDEXING_CHECK_CELL_COMBINATIONS ) need_cell = 1; if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1; - if ( need_cell && (cell == NULL) ) { - ERROR("Altering your DirAx flags because no PDB file was" - " provided.\n"); + if ( need_cell && !cell_has_parameters(cell) ) { + ERROR("Altering your DirAx flags because cell parameters were" + " not provided.\n"); *indm &= ~INDEXING_CHECK_CELL_COMBINATIONS; *indm &= ~INDEXING_CHECK_CELL_AXES; } diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index a93bbfe4..cf0852bf 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -841,9 +841,9 @@ IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1; if ( *indm & INDEXING_USE_LATTICE_TYPE ) need_cell = 1; - if ( need_cell && (cell == NULL) ) { - ERROR("Altering your MOSFLM flags because no PDB file was" - " provided.\n"); + if ( need_cell && !cell_has_parameters(cell) ) { + ERROR("Altering your MOSFLM flags because cell parameters were" + " not provided.\n"); *indm &= ~INDEXING_CHECK_CELL_COMBINATIONS; *indm &= ~INDEXING_CHECK_CELL_AXES; *indm &= ~INDEXING_USE_LATTICE_TYPE; diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c index c29083b5..4a8ef30e 100644 --- a/libcrystfel/src/xds.c +++ b/libcrystfel/src/xds.c @@ -629,9 +629,12 @@ IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell, * but we'd have to decide whether the user just forgot the cell, or * forgot "-nolatt", or whatever. */ if ( ((*indm & INDEXING_USE_LATTICE_TYPE) - || (*indm & INDEXING_USE_CELL_PARAMETERS)) && (cell == NULL) ) { - ERROR("No cell provided. If you wanted to use XDS without " - "prior cell information, use xds-nolatt-nocell.\n"); + || (*indm & INDEXING_USE_CELL_PARAMETERS)) + && !cell_has_parameters(cell) ) + { + ERROR("No cell parameters provided. If you wanted to use XDS " + "without prior cell information, use " + "xds-nolatt-nocell.\n"); return NULL; } diff --git a/src/check_hkl.c b/src/check_hkl.c index f22373b0..e10f5f29 100644 --- a/src/check_hkl.c +++ b/src/check_hkl.c @@ -57,7 +57,7 @@ static void show_help(const char *s) " -h, --help Display this help message.\n" " --version Print CrystFEL version number and exit.\n" " -y, --symmetry=<sym> The symmetry of the input file.\n" -" -p, --pdb=<filename> PDB file to use.\n" +" -p, --pdb=<filename> Unit cell file to use (PDB or CrystFEL format).\n" " --rmin=<res> Low resolution cutoff (1/d in m^-1).\n" " --rmax=<res> High resolution cutoff (1/d in m^-1).\n" " --lowres=<n> Low resolution cutoff in (d in A).\n" @@ -692,7 +692,7 @@ int main(int argc, char *argv[]) RefList *list; Reflection *refl; RefListIterator *iter; - char *pdb = NULL; + char *cellfile = NULL; int rej = 0; float rmin_fix = -1.0; float rmax_fix = -1.0; @@ -751,7 +751,7 @@ int main(int argc, char *argv[]) break; case 'p' : - pdb = strdup(optarg); + cellfile = strdup(optarg); break; case 0 : @@ -833,13 +833,12 @@ int main(int argc, char *argv[]) file = strdup(argv[optind++]); - if ( pdb == NULL ) { - ERROR("You need to provide a PDB file containing" - " the unit cell.\n"); + if ( cellfile == NULL ) { + ERROR("You need to provide a unit cell.\n"); return 1; } - cell = load_cell_from_pdb(pdb); - free(pdb); + cell = load_cell_from_file(cellfile); + free(cellfile); raw_list = read_reflections(file); if ( raw_list == NULL ) { diff --git a/src/compare_hkl.c b/src/compare_hkl.c index e10b93aa..a35e4f3b 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -89,7 +89,7 @@ static void show_help(const char *s) "Compare intensity lists.\n" "\n" " -y, --symmetry=<sym> The symmetry of both the input files.\n" -" -p, --pdb=<filename> PDB file to use.\n" +" -p, --pdb=<filename> Unit cell file to use.\n" " --fom=<FoM> Calculate this figure of merit Choose from:\n" " R1I, R1F, R2, Rsplit, CC, CCstar,\n" " CCano, CRDano, Rano, Rano/Rsplit\n" @@ -916,7 +916,7 @@ int main(int argc, char *argv[]) RefList *list1_raw; RefList *list2_raw; enum fom fom = FOM_R1I; - char *pdb = NULL; + char *cellfile = NULL; float rmin_fix = -1.0; float rmax_fix = -1.0; double rmin, rmax; @@ -974,7 +974,7 @@ int main(int argc, char *argv[]) break; case 'p' : - pdb = strdup(optarg); + cellfile = strdup(optarg); break; case 'u' : @@ -1116,15 +1116,15 @@ int main(int argc, char *argv[]) afile = strdup(argv[optind++]); bfile = strdup(argv[optind]); - if ( pdb == NULL ) { - ERROR("You must provide a PDB file with the unit cell.\n"); + if ( cellfile == NULL ) { + ERROR("You must provide a unit cell.\n"); exit(1); } if ( shell_file == NULL ) shell_file = strdup("shells.dat"); - cell = load_cell_from_pdb(pdb); - free(pdb); + cell = load_cell_from_file(cellfile); + free(cellfile); list1_raw = read_reflections(afile); if ( list1_raw == NULL ) { diff --git a/src/get_hkl.c b/src/get_hkl.c index 3da68ef7..72b73041 100644 --- a/src/get_hkl.c +++ b/src/get_hkl.c @@ -421,7 +421,7 @@ int main(int argc, char *argv[]) char *cutoff_str = NULL; double cutiso = 0.0; float cutn1, cutn2, cutn3; - char *pdb = NULL; + char *cellfile = NULL; char *reindex_str = NULL; SymOpList *reindex = NULL; @@ -487,7 +487,7 @@ int main(int argc, char *argv[]) break; case 'p' : - pdb = strdup(optarg); + cellfile = strdup(optarg); break; case 2 : @@ -715,18 +715,18 @@ int main(int argc, char *argv[]) RefListIterator *iter; UnitCell *cell; - if ( pdb == NULL ) { - ERROR("You must provide a PDB file when using " + if ( cellfile == NULL ) { + ERROR("You must provide a unit cell when using " "--cutoff-angstroms.\n"); return 1; } - cell = load_cell_from_pdb(pdb); + cell = load_cell_from_file(cellfile); if ( cell == NULL ) { - ERROR("Failed to load cell from '%s'\n", pdb); + ERROR("Failed to load cell from '%s'\n", cellfile); return 1; } - free(pdb); + free(cellfile); n = reflist_new(); @@ -762,18 +762,18 @@ int main(int argc, char *argv[]) double csx, csy, csz; double as, bs, cs; - if ( pdb == NULL ) { - ERROR("You must provide a PDB file when using " + if ( cellfile == NULL ) { + ERROR("You must provide a unit cell when using " "--cutoff-angstroms.\n"); return 1; } - cell = load_cell_from_pdb(pdb); + cell = load_cell_from_file(cellfile); if ( cell == NULL ) { - ERROR("Failed to load cell from '%s'\n", pdb); + ERROR("Failed to load cell from '%s'\n", cellfile); return 1; } - free(pdb); + free(cellfile); cell_get_reciprocal(cell, &asx, &asy, &asz, &bsx, &bsy, &bsz, diff --git a/src/indexamajig.c b/src/indexamajig.c index 2ef98c2c..ff4b2243 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -93,8 +93,8 @@ static void show_help(const char *s) " -b, --beam=<file> Get beam parameters from file (provides nominal\n" " wavelength value if no per-shot value is found in\n" " the HDF5 files.\n" -" -p, --pdb=<file> PDB file from which to get the unit cell to match.\n" -" Default: 'molecule.pdb'.\n" +" -p, --pdb=<file> File (PDB or CrystFEL unit cell format) from which\n" +" to get the unit cell. Default: 'molecule.pdb'.\n" " --basename Remove the directory parts of the filenames.\n" " -x, --prefix=<p> Prefix filenames from input file with <p>.\n" " --peaks=<method> Use 'method' for finding peaks. Choose from:\n" @@ -187,7 +187,7 @@ int main(int argc, char *argv[]) IndexingMethod *indm; IndexingPrivate **ipriv; char *indm_str = NULL; - char *pdb = NULL; + char *cellfile = NULL; char *prefix = NULL; char *speaks = NULL; char *toler = NULL; @@ -327,7 +327,7 @@ int main(int argc, char *argv[]) break; case 'p' : - pdb = strdup(optarg); + cellfile = strdup(optarg); break; case 'x' : @@ -606,13 +606,13 @@ int main(int argc, char *argv[]) add_geom_beam_stuff_to_copy_hdf5(iargs.copyme, iargs.det, iargs.beam); - if ( pdb != NULL ) { - iargs.cell = load_cell_from_pdb(pdb); + if ( cellfile != NULL ) { + iargs.cell = load_cell_from_file(cellfile); if ( iargs.cell == NULL ) { - ERROR("Couldn't read unit cell (from %s)\n", pdb); + ERROR("Couldn't read unit cell (from %s)\n", cellfile); return 1; } - free(pdb); + free(cellfile); STATUS("This is what I understood your unit cell to be:\n"); cell_print(iargs.cell); } else { diff --git a/src/partial_sim.c b/src/partial_sim.c index d5933afd..47d8766e 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -624,7 +624,7 @@ int main(int argc, char *argv[]) ERROR("You need to give a PDB file with the unit cell.\n"); return 1; } - cell = load_cell_from_pdb(cellfile); + cell = load_cell_from_file(cellfile); if ( cell == NULL ) { ERROR("Failed to get cell from '%s'\n", cellfile); return 1; diff --git a/src/pattern_sim.c b/src/pattern_sim.c index 0889697d..8e048b4e 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -69,7 +69,7 @@ static void show_help(const char *s) " -h, --help Display this help message.\n" " --version Print CrystFEL version number and exit.\n" "\n" -" -p, --pdb=<file> PDB file from which to get the unit cell.\n" +" -p, --pdb=<file> File from which to get the unit cell.\n" " (The actual Bragg intensities come from the\n" " intensities file)\n" " --gpu Use the GPU to speed up the calculation.\n" @@ -572,7 +572,7 @@ int main(int argc, char *argv[]) free(beamfile); /* Load unit cell */ - input_cell = load_cell_from_pdb(filename); + input_cell = load_cell_from_file(filename); if ( input_cell == NULL ) { exit(1); } diff --git a/src/render_hkl.c b/src/render_hkl.c index bdf1bbe6..1b13be66 100644 --- a/src/render_hkl.c +++ b/src/render_hkl.c @@ -3,11 +3,11 @@ * * Draw pretty renderings of reflection lists * - * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2013 Thomas White <taw@physics.org> + * 2010-2014 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -71,7 +71,7 @@ static void show_help(const char *s) " --zone=<z> Show the <z>th Laue zone.\n" " -o, --output=<filename> Output filename. Default: za.pdf\n" " --boost=<val> Squash colour scale by <val>.\n" -" -p, --pdb=<file> PDB file from which to get the unit cell.\n" +" -p, --pdb=<file> File from which to get the unit cell.\n" " -y, --symmetry=<sym> Expand reflections according to point group <sym>.\n" "\n" " -c, --colscale=<scale> Use the given colour scale. Choose from:\n" @@ -741,7 +741,7 @@ int main(int argc, char *argv[]) int config_sqrt = 0; int config_colkey = 0; int config_zawhinge = 0; - char *pdb = NULL; + char *cellfile = NULL; int r = 0; double boost = 1.0; char *sym_str = NULL; @@ -800,7 +800,7 @@ int main(int argc, char *argv[]) return 0; case 'p' : - pdb = strdup(optarg); + cellfile = strdup(optarg); break; case 'b' : @@ -888,8 +888,8 @@ int main(int argc, char *argv[]) " any longer (I ignored it for you).\n"); } - if ( (pdb == NULL) && !config_colkey ) { - ERROR("You must specify the PDB containing the unit cell.\n"); + if ( (cellfile == NULL) && !config_colkey ) { + ERROR("You must specify the unit cell.\n"); return 1; } @@ -974,9 +974,9 @@ int main(int argc, char *argv[]) infile = argv[optind]; - cell = load_cell_from_pdb(pdb); + cell = load_cell_from_file(cellfile); if ( cell == NULL ) { - ERROR("Couldn't load unit cell from %s\n", pdb); + ERROR("Couldn't load unit cell from %s\n", cellfile); return 1; } list = read_reflections(infile); @@ -993,7 +993,7 @@ int main(int argc, char *argv[]) render_za(cell, list, boost, sym, wght, colscale, rh, rk, rl, dh, dk, dl, outfile, scale_top, zone, &rings); - free(pdb); + free(cellfile); free_symoplist(sym); reflist_free(list); if ( outfile != NULL ) free(outfile); |