aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2014-09-19 16:07:24 +0200
committerThomas White <taw@physics.org>2014-09-19 16:23:39 +0200
commit6a6cb3b4d7f15c234a79ff8421a0ae5c1a1dcb2a (patch)
tree00f6e0da9a8d086af18b0b1f34433bc115c9f206 /libcrystfel
parent2c959daa7a46b99a10dd5a1998b62ccb8def97de (diff)
Introduce CrystFEL unit cell files
Diffstat (limited to 'libcrystfel')
-rw-r--r--libcrystfel/src/cell-utils.c262
-rw-r--r--libcrystfel/src/cell-utils.h9
-rw-r--r--libcrystfel/src/cell.c83
-rw-r--r--libcrystfel/src/cell.h12
-rw-r--r--libcrystfel/src/dirax.c6
-rw-r--r--libcrystfel/src/mosflm.c6
-rw-r--r--libcrystfel/src/xds.c9
7 files changed, 331 insertions, 56 deletions
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;
}