From 6a6cb3b4d7f15c234a79ff8421a0ae5c1a1dcb2a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 19 Sep 2014 16:07:24 +0200 Subject: Introduce CrystFEL unit cell files --- libcrystfel/src/cell-utils.c | 262 ++++++++++++++++++++++++++++++++++++++----- libcrystfel/src/cell-utils.h | 9 +- libcrystfel/src/cell.c | 83 ++++++++++++-- libcrystfel/src/cell.h | 12 +- libcrystfel/src/dirax.c | 6 +- libcrystfel/src/mosflm.c | 6 +- libcrystfel/src/xds.c | 9 +- 7 files changed, 331 insertions(+), 56 deletions(-) (limited to 'libcrystfel') 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 - * 2012 Lorenzo Galli + * 2009-2012,2014 Thomas White + * 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 - * 2012 Lorenzo Galli + * 2009-2013,2014 Thomas White + * 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 - * 2010 Richard Kirian - * 2012 Lorenzo Galli + * 2009-2012,2014 Thomas White + * 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 - * 2010,2012 Richard Kirian - * 2012 Lorenzo Galli + * 2009-2012,2014 Thomas White + * 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; } -- cgit v1.2.3