From ae3a6a04282ffaeaa010edbff1850a522d255054 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 7 Aug 2012 18:20:38 +0200 Subject: Move unit cell utility stuff to separate module --- libcrystfel/Makefile.am | 5 +- libcrystfel/src/cell-utils.c | 863 ++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/cell-utils.h | 58 +++ libcrystfel/src/cell.c | 825 +------------------------------------- libcrystfel/src/cell.h | 22 +- libcrystfel/src/geometry.c | 1 + libcrystfel/src/index.c | 1 + libcrystfel/src/peaks.c | 1 + libcrystfel/src/reax.c | 1 + libcrystfel/src/reflist-utils.c | 1 + 10 files changed, 933 insertions(+), 845 deletions(-) create mode 100644 libcrystfel/src/cell-utils.c create mode 100644 libcrystfel/src/cell-utils.h (limited to 'libcrystfel') diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index d55f6e74..28d483ac 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -7,7 +7,8 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/beam-parameters.c src/geometry.c src/statistics.c \ src/symmetry.c src/stream.c src/peaks.c \ src/reflist-utils.c src/filters.c \ - src/render.c src/index.c src/dirax.c src/mosflm.c + src/render.c src/index.c src/dirax.c src/mosflm.c \ + src/cell-utils.c if HAVE_FFTW libcrystfel_la_SOURCES += src/reax.c @@ -22,7 +23,7 @@ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/hdf5-file.h \ src/geometry.h src/peaks.h src/stream.h \ src/render.h src/index.h src/image.h \ src/filters.h src/dirax.h src/mosflm.h \ - src/index-priv.h src/reax.h + src/index-priv.h src/reax.h src/cell-utils.h INCLUDES = "-I$(top_srcdir)/data" AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c new file mode 100644 index 00000000..c30d9a03 --- /dev/null +++ b/libcrystfel/src/cell-utils.c @@ -0,0 +1,863 @@ +/* + * cell-utils.c + * + * Unit Cell utility functions + * + * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * Copyright © 2012 Lorenzo Galli + * + * Authors: + * 2009-2012 Thomas White + * 2012 Lorenzo Galli + * + * 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 . + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include +#include + +#include "cell.h" +#include "cell-utils.h" +#include "utils.h" +#include "image.h" + + +/* Weighting factor of lengths relative to angles */ +#define LWEIGHT (10.0e-9) + + +/** + * cell_rotate: + * @in: A %UnitCell to rotate + * @quat: A %quaternion + * + * Rotate a %UnitCell using a %quaternion. + * + * Returns: a newly allocated rotated copy of @in. + * + */ +UnitCell *cell_rotate(UnitCell *in, struct quaternion quat) +{ + struct rvec a, b, c; + struct rvec an, bn, cn; + UnitCell *out = cell_new_from_cell(in); + + cell_get_cartesian(in, &a.u, &a.v, &a.w, + &b.u, &b.v, &b.w, + &c.u, &c.v, &c.w); + + an = quat_rot(a, quat); + bn = quat_rot(b, quat); + cn = quat_rot(c, quat); + + cell_set_cartesian(out, an.u, an.v, an.w, + bn.u, bn.v, bn.w, + cn.u, cn.v, cn.w); + + return out; +} + + +static const char *str_lattice(LatticeType l) +{ + switch ( l ) + { + case L_TRICLINIC : return "triclinic"; + case L_MONOCLINIC : return "monoclinic"; + case L_ORTHORHOMBIC : return "orthorhombic"; + case L_TETRAGONAL : return "tetragonal"; + case L_RHOMBOHEDRAL : return "rhombohedral"; + case L_HEXAGONAL : return "hexagonal"; + case L_CUBIC : return "cubic"; + } + + return "unknown lattice"; +} + + +void cell_print(UnitCell *cell) +{ + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double a, b, c, alpha, beta, gamma; + double ax, ay, az, bx, by, bz, cx, cy, cz; + + STATUS("%s %c\n", str_lattice(cell_get_lattice_type(cell)), + cell_get_centering(cell)); + + cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); + + STATUS(" a b c alpha beta gamma\n"); + STATUS("%5.2f %5.2f %5.2f nm %6.2f %6.2f %6.2f deg\n", + a*1e9, b*1e9, c*1e9, + rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); + + 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); + + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &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("Point group: %s\n", cell_get_pointgroup(cell)); + STATUS("Cell representation is %s.\n", cell_rep(cell)); +} + + +#define MAX_CAND (1024) + +static int right_handed(struct rvec a, struct rvec b, struct rvec c) +{ + struct rvec aCb; + double aCb_dot_c; + + /* "a" cross "b" */ + aCb.u = a.v*b.w - a.w*b.v; + aCb.v = - (a.u*b.w - a.w*b.u); + aCb.w = a.u*b.v - a.v*b.u; + + /* "a cross b" dot "c" */ + aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w; + + if ( aCb_dot_c > 0.0 ) return 1; + return 0; +} + + +struct cvec { + struct rvec vec; + float na; + float nb; + float nc; + float fom; +}; + + +static int same_vector(struct cvec a, struct cvec b) +{ + if ( a.na != b.na ) return 0; + if ( a.nb != b.nb ) return 0; + if ( a.nc != b.nc ) return 0; + return 1; +} + + +/* Attempt to make 'cell' fit into 'template' somehow */ +UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, + const float *tols, int reduce) +{ + signed int n1l, n2l, n3l; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + int i, j; + double lengths[3]; + double angles[3]; + struct cvec *cand[3]; + UnitCell *new_cell = NULL; + float best_fom = +999999999.9; /* Large number.. */ + int ncand[3] = {0,0,0}; + signed int ilow, ihigh; + float angtol = deg2rad(tols[3]); + + if ( cell_get_reciprocal(template, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz) ) { + ERROR("Couldn't get reciprocal cell for template.\n"); + return NULL; + } + + lengths[0] = modulus(asx, asy, asz); + lengths[1] = modulus(bsx, bsy, bsz); + lengths[2] = modulus(csx, csy, csz); + + angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz); + angles[1] = angle_between(asx, asy, asz, csx, csy, csz); + angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz); + + cand[0] = malloc(MAX_CAND*sizeof(struct cvec)); + cand[1] = malloc(MAX_CAND*sizeof(struct cvec)); + cand[2] = malloc(MAX_CAND*sizeof(struct cvec)); + + if ( cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz) ) { + ERROR("Couldn't get reciprocal cell.\n"); + return NULL; + } + + if ( reduce ) { + ilow = -2; ihigh = 4; + } else { + ilow = 0; ihigh = 1; + } + + /* Negative values mean 1/n, positive means n, zero means zero */ + for ( n1l=ilow; n1l<=ihigh; n1l++ ) { + for ( n2l=ilow; n2l<=ihigh; n2l++ ) { + for ( n3l=ilow; n3l<=ihigh; n3l++ ) { + + float n1, n2, n3; + signed int b1, b2, b3; + + n1 = (n1l>=0) ? (n1l) : (1.0/n1l); + n2 = (n2l>=0) ? (n2l) : (1.0/n2l); + n3 = (n3l>=0) ? (n3l) : (1.0/n3l); + + if ( !reduce ) { + if ( n1l + n2l + n3l > 1 ) continue; + } + + /* 'bit' values can be +1 or -1 */ + for ( b1=-1; b1<=1; b1+=2 ) { + for ( b2=-1; b2<=1; b2+=2 ) { + for ( b3=-1; b3<=1; b3+=2 ) { + + double tx, ty, tz; + double tlen; + int i; + + n1 *= b1; n2 *= b2; n3 *= b3; + + tx = n1*asx + n2*bsx + n3*csx; + ty = n1*asy + n2*bsy + n3*csy; + tz = n1*asz + n2*bsz + n3*csz; + tlen = modulus(tx, ty, tz); + + /* Test modulus for agreement with moduli of template */ + for ( i=0; i<3; i++ ) { + + if ( !within_tolerance(lengths[i], tlen, + tols[i]) ) + { + continue; + } + + if ( ncand[i] == MAX_CAND ) { + ERROR("Too many cell candidates - "); + ERROR("consider tightening the unit "); + ERROR("cell tolerances.\n"); + } else { + + double fom; + + fom = fabs(lengths[i] - tlen); + + cand[i][ncand[i]].vec.u = tx; + cand[i][ncand[i]].vec.v = ty; + cand[i][ncand[i]].vec.w = tz; + cand[i][ncand[i]].na = n1; + cand[i][ncand[i]].nb = n2; + cand[i][ncand[i]].nc = n3; + cand[i][ncand[i]].fom = fom; + + ncand[i]++; + + } + } + + } + } + } + + } + } + } + + if ( verbose ) { + STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]); + } + + for ( i=0; i angtol ) continue; + + fom1 = fabs(ang - angles[2]); + + for ( k=0; k angtol ) continue; + + fom2 = fom1 + fabs(ang - angles[1]); + + /* Finally, the angle between the current candidate for + * axis 1 and the kth candidate for axis 2 */ + ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v, + cand[1][j].vec.w, cand[2][k].vec.u, + cand[2][k].vec.v, cand[2][k].vec.w); + + /* ... it should be angle 0 ... */ + if ( fabs(ang - angles[0]) > angtol ) continue; + + /* Unit cell must be right-handed */ + if ( !right_handed(cand[0][i].vec, cand[1][j].vec, + cand[2][k].vec) ) continue; + + fom3 = fom2 + fabs(ang - angles[0]); + fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom + + cand[2][k].fom); + + if ( fom3 < best_fom ) { + if ( new_cell != NULL ) free(new_cell); + new_cell = cell_new_from_reciprocal_axes( + cand[0][i].vec, cand[1][j].vec, + cand[2][k].vec); + best_fom = fom3; + } + + } + + } + } + + free(cand[0]); + free(cand[1]); + free(cand[2]); + + return new_cell; +} + + + +UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template) +{ + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + int i; + double lengths[3]; + int used[3]; + struct rvec real_a, real_b, real_c; + struct rvec params[3]; + double alen, blen; + float ltl = 5.0; /* percent */ + int have_real_a; + int have_real_b; + int have_real_c; + + /* Get the lengths to match */ + if ( cell_get_cartesian(template, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz) ) + { + ERROR("Couldn't get cell for template.\n"); + return NULL; + } + alen = modulus(ax, ay, az); + blen = modulus(bx, by, bz); + + /* Get the lengths from the cell and turn them into anonymous vectors */ + if ( cell_get_cartesian(cell, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz) ) + { + ERROR("Couldn't get cell.\n"); + return NULL; + } + lengths[0] = modulus(ax, ay, az); + lengths[1] = modulus(bx, by, bz); + lengths[2] = modulus(cx, cy, cz); + used[0] = 0; used[1] = 0; used[2] = 0; + params[0].u = ax; params[0].v = ay; params[0].w = az; + params[1].u = bx; params[1].v = by; params[1].w = bz; + params[2].u = cx; params[2].v = cy; params[2].w = cz; + + real_a.u = 0.0; real_a.v = 0.0; real_a.w = 0.0; + real_b.u = 0.0; real_b.v = 0.0; real_b.w = 0.0; + real_c.u = 0.0; real_c.v = 0.0; real_c.w = 0.0; + + /* Check each vector against a and b */ + have_real_a = 0; + have_real_b = 0; + for ( i=0; i<3; i++ ) { + if ( within_tolerance(lengths[i], alen, ltl) + && !used[i] && !have_real_a ) + { + used[i] = 1; + memcpy(&real_a, ¶ms[i], sizeof(struct rvec)); + have_real_a = 1; + } + if ( within_tolerance(lengths[i], blen, ltl) + && !used[i] && !have_real_b ) + { + used[i] = 1; + memcpy(&real_b, ¶ms[i], sizeof(struct rvec)); + have_real_b = 1; + } + } + + /* Have we matched both a and b? */ + if ( !(have_real_a && have_real_b) ) return NULL; + + /* "c" is "the other one" */ + have_real_c = 0; + for ( i=0; i<3; i++ ) { + if ( !used[i] ) { + memcpy(&real_c, ¶ms[i], sizeof(struct rvec)); + have_real_c = 1; + } + } + + if ( !have_real_c ) { + ERROR("Huh? Couldn't find the third vector.\n"); + ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]); + return NULL; + } + + /* Flip c if not right-handed */ + if ( !right_handed(real_a, real_b, real_c) ) { + real_c.u = -real_c.u; + real_c.v = -real_c.v; + real_c.w = -real_c.w; + } + + return cell_new_from_direct_axes(real_a, real_b, real_c); +} + + +/* Return sin(theta)/lambda = 1/2d. Multiply by two if you want 1/d */ +double resolution(UnitCell *cell, signed int h, signed int k, signed int l) +{ + double a, b, c, alpha, beta, gamma; + + cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); + + const double Vsq = a*a*b*b*c*c*(1 - cos(alpha)*cos(alpha) + - cos(beta)*cos(beta) + - cos(gamma)*cos(gamma) + + 2*cos(alpha)*cos(beta)*cos(gamma) ); + + const double S11 = b*b*c*c*sin(alpha)*sin(alpha); + const double S22 = a*a*c*c*sin(beta)*sin(beta); + const double S33 = a*a*b*b*sin(gamma)*sin(gamma); + const double S12 = a*b*c*c*(cos(alpha)*cos(beta) - cos(gamma)); + const double S23 = a*a*b*c*(cos(beta)*cos(gamma) - cos(alpha)); + const double S13 = a*b*b*c*(cos(gamma)*cos(alpha) - cos(beta)); + + const double brackets = S11*h*h + S22*k*k + S33*l*l + + 2*S12*h*k + 2*S23*k*l + 2*S13*h*l; + const double oneoverdsq = brackets / Vsq; + const double oneoverd = sqrt(oneoverdsq); + + return oneoverd / 2; +} + + +static void determine_lattice(UnitCell *cell, + const char *as, const char *bs, const char *cs, + const char *als, const char *bes, const char *gas) +{ + int n_right; + + /* Rhombohedral or cubic? */ + if ( (strcmp(as, bs) == 0) && (strcmp(as, cs) == 0) ) { + + if ( (strcmp(als, " 90.00") == 0) + && (strcmp(bes, " 90.00") == 0) + && (strcmp(gas, " 90.00") == 0) ) + { + /* Cubic. Unique axis irrelevant. */ + cell_set_lattice_type(cell, L_CUBIC); + return; + } + + if ( (strcmp(als, bes) == 0) && (strcmp(als, gas) == 0) ) { + /* Rhombohedral. Unique axis irrelevant. */ + cell_set_lattice_type(cell, L_RHOMBOHEDRAL); + return; + } + + } + + if ( (strcmp(als, " 90.00") == 0) + && (strcmp(bes, " 90.00") == 0) + && (strcmp(gas, " 90.00") == 0) ) + { + if ( strcmp(bs, cs) == 0 ) { + /* Tetragonal, unique axis a */ + cell_set_lattice_type(cell, L_TETRAGONAL); + cell_set_unique_axis(cell, 'a'); + return; + } + + if ( strcmp(as, cs) == 0 ) { + /* Tetragonal, unique axis b */ + cell_set_lattice_type(cell, L_TETRAGONAL); + cell_set_unique_axis(cell, 'b'); + return; + } + + if ( strcmp(as, bs) == 0 ) { + /* Tetragonal, unique axis c */ + cell_set_lattice_type(cell, L_TETRAGONAL); + cell_set_unique_axis(cell, 'c'); + return; + } + + /* Orthorhombic. Unique axis irrelevant, but point group + * can have different orientations. */ + cell_set_lattice_type(cell, L_ORTHORHOMBIC); + return; + } + + n_right = 0; + if ( strcmp(als, " 90.00") == 0 ) n_right++; + if ( strcmp(bes, " 90.00") == 0 ) n_right++; + if ( strcmp(gas, " 90.00") == 0 ) n_right++; + + /* Hexgonal or monoclinic? */ + if ( n_right == 2 ) { + + if ( (strcmp(als, " 120.00") == 0) + && (strcmp(bs, cs) == 0) ) + { + /* Hexagonal, unique axis a */ + cell_set_lattice_type(cell, L_HEXAGONAL); + cell_set_unique_axis(cell, 'a'); + return; + } + + if ( (strcmp(bes, " 120.00") == 0) + && (strcmp(as, cs) == 0) ) + { + /* Hexagonal, unique axis b */ + cell_set_lattice_type(cell, L_HEXAGONAL); + cell_set_unique_axis(cell, 'b'); + return; + } + + if ( (strcmp(gas, " 120.00") == 0) + && (strcmp(as, bs) == 0) ) + { + /* Hexagonal, unique axis c */ + cell_set_lattice_type(cell, L_HEXAGONAL); + cell_set_unique_axis(cell, 'c'); + return; + } + + if ( strcmp(als, " 90.00") != 0 ) { + /* Monoclinic, unique axis a */ + cell_set_lattice_type(cell, L_MONOCLINIC); + cell_set_unique_axis(cell, 'a'); + return; + } + + if ( strcmp(bes, " 90.00") != 0 ) { + /* Monoclinic, unique axis b */ + cell_set_lattice_type(cell, L_MONOCLINIC); + cell_set_unique_axis(cell, 'b'); + return; + } + + if ( strcmp(gas, " 90.00") != 0 ) { + /* Monoclinic, unique axis c */ + cell_set_lattice_type(cell, L_MONOCLINIC); + cell_set_unique_axis(cell, 'c'); + return; + } + } + + /* Triclinic, unique axis irrelevant. */ + cell_set_lattice_type(cell, L_TRICLINIC); +} + + +UnitCell *load_cell_from_pdb(const char *filename) +{ + FILE *fh; + char *rval; + UnitCell *cell = NULL; + + fh = fopen(filename, "r"); + if ( fh == NULL ) { + ERROR("Couldn't open '%s'\n", filename); + return NULL; + } + + do { + + char line[1024]; + + rval = fgets(line, 1023, fh); + + if ( strncmp(line, "CRYST1", 6) == 0 ) { + + float a, b, c, al, be, ga; + char as[10], bs[10], cs[10]; + char als[8], bes[8], gas[8]; + int r; + + memcpy(as, line+6, 9); as[9] = '\0'; + memcpy(bs, line+15, 9); bs[9] = '\0'; + memcpy(cs, line+24, 9); cs[9] = '\0'; + memcpy(als, line+33, 7); als[7] = '\0'; + memcpy(bes, line+40, 7); bes[7] = '\0'; + memcpy(gas, line+47, 7); gas[7] = '\0'; + + STATUS("'%s' '%s' '%s'\n", as, bs, cs); + STATUS("'%s' '%s' '%s'\n", als, bes, gas); + + r = sscanf(as, "%f", &a); + r += sscanf(bs, "%f", &b); + r += sscanf(cs, "%f", &c); + r += sscanf(als, "%f", &al); + r += sscanf(bes, "%f", &be); + r += sscanf(gas, "%f", &ga); + + if ( r != 6 ) { + STATUS("Couldn't understand CRYST1 line.\n"); + continue; + } + + cell = cell_new_from_parameters(a*1e-10, + b*1e-10, c*1e-10, + deg2rad(al), + deg2rad(be), + deg2rad(ga)); + + determine_lattice(cell, as, bs, cs, als, bes, gas); + + if ( strlen(line) > 65 ) { + cell_set_centering(cell, line[55]); + } else { + cell_set_pointgroup(cell, "1"); + ERROR("CRYST1 line without centering.\n"); + } + + break; /* Done */ + } + + } while ( rval != NULL ); + + fclose(fh); + + /* FIXME: Turn "H" centered cells into "R" cells */ + + validate_cell(cell); + + return cell; +} + + +/* Force the linker to bring in CBLAS to make GSL happy */ +void cell_fudge_gslcblas() +{ + STATUS("%p\n", cblas_sgemm); +} + + +UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot) +{ + UnitCell *out; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double xnew, ynew, znew; + + cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy, + &bsz, &csx, &csy, &csz); + + /* Rotate by "omega" about +z (parallel to c* and c unless triclinic) */ + xnew = asx*cos(omega) + asy*sin(omega); + ynew = -asx*sin(omega) + asy*cos(omega); + znew = asz; + asx = xnew; asy = ynew; asz = znew; + xnew = bsx*cos(omega) + bsy*sin(omega); + ynew = -bsx*sin(omega) + bsy*cos(omega); + znew = bsz; + bsx = xnew; bsy = ynew; bsz = znew; + xnew = csx*cos(omega) + csy*sin(omega); + ynew = -csx*sin(omega) + csy*cos(omega); + znew = csz; + csx = xnew; csy = ynew; csz = znew; + + /* Rotate by "phi" about +x (not parallel to anything specific) */ + xnew = asx; + ynew = asy*cos(phi) + asz*sin(phi); + znew = -asy*sin(phi) + asz*cos(phi); + asx = xnew; asy = ynew; asz = znew; + xnew = bsx; + ynew = bsy*cos(phi) + bsz*sin(phi); + znew = -bsy*sin(phi) + bsz*cos(phi); + bsx = xnew; bsy = ynew; bsz = znew; + xnew = csx; + ynew = csy*cos(phi) + csz*sin(phi); + znew = -csy*sin(phi) + csz*cos(phi); + csx = xnew; csy = ynew; csz = znew; + + /* Rotate by "rot" about the new +z (in-plane rotation) */ + xnew = asx*cos(rot) + asy*sin(rot); + ynew = -asx*sin(rot) + asy*cos(rot); + znew = asz; + asx = xnew; asy = ynew; asz = znew; + xnew = bsx*cos(rot) + bsy*sin(rot); + ynew = -bsx*sin(rot) + bsy*cos(rot); + znew = bsz; + bsx = xnew; bsy = ynew; bsz = znew; + xnew = csx*cos(rot) + csy*sin(rot); + ynew = -csx*sin(rot) + csy*cos(rot); + znew = csz; + csx = xnew; csy = ynew; csz = znew; + + out = cell_new_from_cell(in); + cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + + return out; +} + + +int cell_is_sensible(UnitCell *cell) +{ + double a, b, c, al, be, ga; + + cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); + if ( al + be + ga >= 2.0*M_PI ) return 0; + if ( al + be - ga >= 2.0*M_PI ) return 0; + if ( al - be + ga >= 2.0*M_PI ) return 0; + if ( - al + be + ga >= 2.0*M_PI ) return 0; + if ( al + be + ga <= 0.0 ) return 0; + if ( al + be - ga <= 0.0 ) return 0; + if ( al - be + ga <= 0.0 ) return 0; + if ( - al + be + ga <= 0.0 ) return 0; + if ( isnan(al) ) return 0; + if ( isnan(be) ) return 0; + if ( isnan(ga) ) return 0; + return 1; +} + + +static int bravais_lattice(UnitCell *cell) +{ + LatticeType lattice = cell_get_lattice_type(cell); + char centering = cell_get_centering(cell); + char ua = cell_get_unique_axis(cell); + + switch ( centering ) + { + case 'P' : + return 1; + + case 'A' : + case 'B' : + case 'C' : + if ( (lattice != L_MONOCLINIC) + && (lattice != L_ORTHORHOMBIC) ) + { + return 0; + } + if ( (ua=='a') && (centering=='A') ) return 1; + if ( (ua=='b') && (centering=='B') ) return 1; + if ( (ua=='c') && (centering=='C') ) return 1; + return 0; + + case 'I' : + if ( (lattice == L_ORTHORHOMBIC) + || (lattice == L_TETRAGONAL) + || (lattice == L_CUBIC) ) + { + return 1; + } + return 0; + + case 'F' : + if ( (lattice == L_ORTHORHOMBIC) || (lattice == L_CUBIC) ) { + return 1; + } + return 0; + + case 'H' : + if ( lattice == L_HEXAGONAL ) return 1; + return 0; + + default : + return 0; + } +} + + +/** + * validate_cell: + * @cell: A %UnitCell to validate + * + * Perform some checks for crystallographic validity @cell, such as that the + * lattice is a conventional Bravais lattice. + * Warnings are printied if any of the checks are failed. + * + */ +void validate_cell(UnitCell *cell) +{ + int err = 0; + + if ( !cell_is_sensible(cell) ) { + ERROR("Warning: Unit cell parameters are not sensible.\n"); + err = 1; + } + + if ( !bravais_lattice(cell) ) { + ERROR("Warning: Unit cell is not a conventional Bravais" + " lattice.\n"); + err = 1; + } + + if ( err ) cell_print(cell); +} diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h new file mode 100644 index 00000000..8d8cfc6d --- /dev/null +++ b/libcrystfel/src/cell-utils.h @@ -0,0 +1,58 @@ +/* + * cell-utils.h + * + * Unit Cell utility functions + * + * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * Copyright © 2012 Lorenzo Galli + * + * Authors: + * 2009-2012 Thomas White + * 2012 Lorenzo Galli + * + * 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 . + * + */ + +#ifndef CELL_UTILS_H +#define CELL_UTILS_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +extern double resolution(UnitCell *cell, + signed int h, signed int k, signed int l); + +extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat); +extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi, + double rot); + +extern void cell_print(UnitCell *cell); + +extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose, + const float *ltl, int reduce); + +extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *tempcell); + +extern UnitCell *load_cell_from_pdb(const char *filename); + +extern int cell_is_sensible(UnitCell *cell); + +extern void validate_cell(UnitCell *cell); + +#endif /* CELL_UTILS_H */ diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 8f380128..64bcd3f4 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -1,7 +1,7 @@ /* * cell.c * - * Unit Cell Calculations + * A class representing a unit cell * * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. @@ -60,10 +60,6 @@ */ - -/* Weighting factor of lengths relative to angles */ -#define LWEIGHT (10.0e-9) - typedef enum { CELL_REP_CRYST, CELL_REP_CART, @@ -603,11 +599,7 @@ char cell_get_unique_axis(UnitCell *cell) } - - -/********************************* Utilities **********************************/ - -static const char *cell_rep(UnitCell *cell) +const char *cell_rep(UnitCell *cell) { switch ( cell->rep ) { @@ -624,816 +616,3 @@ static const char *cell_rep(UnitCell *cell) return "unknown"; } - - -/** - * cell_rotate: - * @in: A %UnitCell to rotate - * @quat: A %quaternion - * - * Rotate a %UnitCell using a %quaternion. - * - * Returns: a newly allocated rotated copy of @in. - * - */ -UnitCell *cell_rotate(UnitCell *in, struct quaternion quat) -{ - struct rvec a, b, c; - struct rvec an, bn, cn; - UnitCell *out = cell_new_from_cell(in); - - cell_get_cartesian(in, &a.u, &a.v, &a.w, - &b.u, &b.v, &b.w, - &c.u, &c.v, &c.w); - - an = quat_rot(a, quat); - bn = quat_rot(b, quat); - cn = quat_rot(c, quat); - - cell_set_cartesian(out, an.u, an.v, an.w, - bn.u, bn.v, bn.w, - cn.u, cn.v, cn.w); - - return out; -} - - -static const char *str_lattice(LatticeType l) -{ - switch ( l ) - { - case L_TRICLINIC : return "triclinic"; - case L_MONOCLINIC : return "monoclinic"; - case L_ORTHORHOMBIC : return "orthorhombic"; - case L_TETRAGONAL : return "tetragonal"; - case L_RHOMBOHEDRAL : return "rhombohedral"; - case L_HEXAGONAL : return "hexagonal"; - case L_CUBIC : return "cubic"; - } - - return "unknown lattice"; -} - - -void cell_print(UnitCell *cell) -{ - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double a, b, c, alpha, beta, gamma; - double ax, ay, az, bx, by, bz, cx, cy, cz; - - STATUS("%s %c\n", str_lattice(cell_get_lattice_type(cell)), - cell_get_centering(cell)); - - cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); - - STATUS(" a b c alpha beta gamma\n"); - STATUS("%5.2f %5.2f %5.2f nm %6.2f %6.2f %6.2f deg\n", - a*1e9, b*1e9, c*1e9, - rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); - - 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); - - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &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("Point group: %s\n", cell_get_pointgroup(cell)); - STATUS("Cell representation is %s.\n", cell_rep(cell)); -} - - -#define MAX_CAND (1024) - -static int right_handed(struct rvec a, struct rvec b, struct rvec c) -{ - struct rvec aCb; - double aCb_dot_c; - - /* "a" cross "b" */ - aCb.u = a.v*b.w - a.w*b.v; - aCb.v = - (a.u*b.w - a.w*b.u); - aCb.w = a.u*b.v - a.v*b.u; - - /* "a cross b" dot "c" */ - aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w; - - if ( aCb_dot_c > 0.0 ) return 1; - return 0; -} - - -struct cvec { - struct rvec vec; - float na; - float nb; - float nc; - float fom; -}; - - -static int same_vector(struct cvec a, struct cvec b) -{ - if ( a.na != b.na ) return 0; - if ( a.nb != b.nb ) return 0; - if ( a.nc != b.nc ) return 0; - return 1; -} - - -/* Attempt to make 'cell' fit into 'template' somehow */ -UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, - const float *tols, int reduce) -{ - signed int n1l, n2l, n3l; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - int i, j; - double lengths[3]; - double angles[3]; - struct cvec *cand[3]; - UnitCell *new_cell = NULL; - float best_fom = +999999999.9; /* Large number.. */ - int ncand[3] = {0,0,0}; - signed int ilow, ihigh; - float angtol = deg2rad(tols[3]); - - if ( cell_get_reciprocal(template, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz) ) { - ERROR("Couldn't get reciprocal cell for template.\n"); - return NULL; - } - - lengths[0] = modulus(asx, asy, asz); - lengths[1] = modulus(bsx, bsy, bsz); - lengths[2] = modulus(csx, csy, csz); - - angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz); - angles[1] = angle_between(asx, asy, asz, csx, csy, csz); - angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz); - - cand[0] = malloc(MAX_CAND*sizeof(struct cvec)); - cand[1] = malloc(MAX_CAND*sizeof(struct cvec)); - cand[2] = malloc(MAX_CAND*sizeof(struct cvec)); - - if ( cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz) ) { - ERROR("Couldn't get reciprocal cell.\n"); - return NULL; - } - - if ( reduce ) { - ilow = -2; ihigh = 4; - } else { - ilow = 0; ihigh = 1; - } - - /* Negative values mean 1/n, positive means n, zero means zero */ - for ( n1l=ilow; n1l<=ihigh; n1l++ ) { - for ( n2l=ilow; n2l<=ihigh; n2l++ ) { - for ( n3l=ilow; n3l<=ihigh; n3l++ ) { - - float n1, n2, n3; - signed int b1, b2, b3; - - n1 = (n1l>=0) ? (n1l) : (1.0/n1l); - n2 = (n2l>=0) ? (n2l) : (1.0/n2l); - n3 = (n3l>=0) ? (n3l) : (1.0/n3l); - - if ( !reduce ) { - if ( n1l + n2l + n3l > 1 ) continue; - } - - /* 'bit' values can be +1 or -1 */ - for ( b1=-1; b1<=1; b1+=2 ) { - for ( b2=-1; b2<=1; b2+=2 ) { - for ( b3=-1; b3<=1; b3+=2 ) { - - double tx, ty, tz; - double tlen; - int i; - - n1 *= b1; n2 *= b2; n3 *= b3; - - tx = n1*asx + n2*bsx + n3*csx; - ty = n1*asy + n2*bsy + n3*csy; - tz = n1*asz + n2*bsz + n3*csz; - tlen = modulus(tx, ty, tz); - - /* Test modulus for agreement with moduli of template */ - for ( i=0; i<3; i++ ) { - - if ( !within_tolerance(lengths[i], tlen, - tols[i]) ) - { - continue; - } - - if ( ncand[i] == MAX_CAND ) { - ERROR("Too many cell candidates - "); - ERROR("consider tightening the unit "); - ERROR("cell tolerances.\n"); - } else { - - double fom; - - fom = fabs(lengths[i] - tlen); - - cand[i][ncand[i]].vec.u = tx; - cand[i][ncand[i]].vec.v = ty; - cand[i][ncand[i]].vec.w = tz; - cand[i][ncand[i]].na = n1; - cand[i][ncand[i]].nb = n2; - cand[i][ncand[i]].nc = n3; - cand[i][ncand[i]].fom = fom; - - ncand[i]++; - - } - - } - - } - } - } - - } - } - } - - if ( verbose ) { - STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]); - } - - for ( i=0; i angtol ) continue; - - fom1 = fabs(ang - angles[2]); - - for ( k=0; k angtol ) continue; - - fom2 = fom1 + fabs(ang - angles[1]); - - /* Finally, the angle between the current candidate for - * axis 1 and the kth candidate for axis 2 */ - ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v, - cand[1][j].vec.w, cand[2][k].vec.u, - cand[2][k].vec.v, cand[2][k].vec.w); - - /* ... it should be angle 0 ... */ - if ( fabs(ang - angles[0]) > angtol ) continue; - - /* Unit cell must be right-handed */ - if ( !right_handed(cand[0][i].vec, cand[1][j].vec, - cand[2][k].vec) ) continue; - - fom3 = fom2 + fabs(ang - angles[0]); - fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom - + cand[2][k].fom); - - if ( fom3 < best_fom ) { - if ( new_cell != NULL ) free(new_cell); - new_cell = cell_new_from_reciprocal_axes( - cand[0][i].vec, cand[1][j].vec, - cand[2][k].vec); - best_fom = fom3; - } - - } - - } - } - - free(cand[0]); - free(cand[1]); - free(cand[2]); - - return new_cell; -} - - - -UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template) -{ - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - int i; - double lengths[3]; - int used[3]; - struct rvec real_a, real_b, real_c; - struct rvec params[3]; - double alen, blen; - float ltl = 5.0; /* percent */ - int have_real_a; - int have_real_b; - int have_real_c; - - /* Get the lengths to match */ - if ( cell_get_cartesian(template, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz) ) - { - ERROR("Couldn't get cell for template.\n"); - return NULL; - } - alen = modulus(ax, ay, az); - blen = modulus(bx, by, bz); - - /* Get the lengths from the cell and turn them into anonymous vectors */ - if ( cell_get_cartesian(cell, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz) ) - { - ERROR("Couldn't get cell.\n"); - return NULL; - } - lengths[0] = modulus(ax, ay, az); - lengths[1] = modulus(bx, by, bz); - lengths[2] = modulus(cx, cy, cz); - used[0] = 0; used[1] = 0; used[2] = 0; - params[0].u = ax; params[0].v = ay; params[0].w = az; - params[1].u = bx; params[1].v = by; params[1].w = bz; - params[2].u = cx; params[2].v = cy; params[2].w = cz; - - real_a.u = 0.0; real_a.v = 0.0; real_a.w = 0.0; - real_b.u = 0.0; real_b.v = 0.0; real_b.w = 0.0; - real_c.u = 0.0; real_c.v = 0.0; real_c.w = 0.0; - - /* Check each vector against a and b */ - have_real_a = 0; - have_real_b = 0; - for ( i=0; i<3; i++ ) { - if ( within_tolerance(lengths[i], alen, ltl) - && !used[i] && !have_real_a ) - { - used[i] = 1; - memcpy(&real_a, ¶ms[i], sizeof(struct rvec)); - have_real_a = 1; - } - if ( within_tolerance(lengths[i], blen, ltl) - && !used[i] && !have_real_b ) - { - used[i] = 1; - memcpy(&real_b, ¶ms[i], sizeof(struct rvec)); - have_real_b = 1; - } - } - - /* Have we matched both a and b? */ - if ( !(have_real_a && have_real_b) ) return NULL; - - /* "c" is "the other one" */ - have_real_c = 0; - for ( i=0; i<3; i++ ) { - if ( !used[i] ) { - memcpy(&real_c, ¶ms[i], sizeof(struct rvec)); - have_real_c = 1; - } - } - - if ( !have_real_c ) { - ERROR("Huh? Couldn't find the third vector.\n"); - ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]); - return NULL; - } - - /* Flip c if not right-handed */ - if ( !right_handed(real_a, real_b, real_c) ) { - real_c.u = -real_c.u; - real_c.v = -real_c.v; - real_c.w = -real_c.w; - } - - return cell_new_from_direct_axes(real_a, real_b, real_c); -} - - -/* Return sin(theta)/lambda = 1/2d. Multiply by two if you want 1/d */ -double resolution(UnitCell *cell, signed int h, signed int k, signed int l) -{ - double a, b, c, alpha, beta, gamma; - - cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); - - const double Vsq = a*a*b*b*c*c*(1 - cos(alpha)*cos(alpha) - - cos(beta)*cos(beta) - - cos(gamma)*cos(gamma) - + 2*cos(alpha)*cos(beta)*cos(gamma) ); - - const double S11 = b*b*c*c*sin(alpha)*sin(alpha); - const double S22 = a*a*c*c*sin(beta)*sin(beta); - const double S33 = a*a*b*b*sin(gamma)*sin(gamma); - const double S12 = a*b*c*c*(cos(alpha)*cos(beta) - cos(gamma)); - const double S23 = a*a*b*c*(cos(beta)*cos(gamma) - cos(alpha)); - const double S13 = a*b*b*c*(cos(gamma)*cos(alpha) - cos(beta)); - - const double brackets = S11*h*h + S22*k*k + S33*l*l - + 2*S12*h*k + 2*S23*k*l + 2*S13*h*l; - const double oneoverdsq = brackets / Vsq; - const double oneoverd = sqrt(oneoverdsq); - - return oneoverd / 2; -} - - -static void determine_lattice(UnitCell *cell, - const char *as, const char *bs, const char *cs, - const char *als, const char *bes, const char *gas) -{ - int n_right; - - /* Rhombohedral or cubic? */ - if ( (strcmp(as, bs) == 0) && (strcmp(as, cs) == 0) ) { - - if ( (strcmp(als, " 90.00") == 0) - && (strcmp(bes, " 90.00") == 0) - && (strcmp(gas, " 90.00") == 0) ) - { - /* Cubic. Unique axis irrelevant. */ - cell_set_lattice_type(cell, L_CUBIC); - return; - } - - if ( (strcmp(als, bes) == 0) && (strcmp(als, gas) == 0) ) { - /* Rhombohedral. Unique axis irrelevant. */ - cell_set_lattice_type(cell, L_RHOMBOHEDRAL); - return; - } - - } - - if ( (strcmp(als, " 90.00") == 0) - && (strcmp(bes, " 90.00") == 0) - && (strcmp(gas, " 90.00") == 0) ) - { - if ( strcmp(bs, cs) == 0 ) { - /* Tetragonal, unique axis a */ - cell_set_lattice_type(cell, L_TETRAGONAL); - cell_set_unique_axis(cell, 'a'); - return; - } - - if ( strcmp(as, cs) == 0 ) { - /* Tetragonal, unique axis b */ - cell_set_lattice_type(cell, L_TETRAGONAL); - cell_set_unique_axis(cell, 'b'); - return; - } - - if ( strcmp(as, bs) == 0 ) { - /* Tetragonal, unique axis c */ - cell_set_lattice_type(cell, L_TETRAGONAL); - cell_set_unique_axis(cell, 'c'); - return; - } - - /* Orthorhombic. Unique axis irrelevant, but point group - * can have different orientations. */ - cell_set_lattice_type(cell, L_ORTHORHOMBIC); - return; - } - - n_right = 0; - if ( strcmp(als, " 90.00") == 0 ) n_right++; - if ( strcmp(bes, " 90.00") == 0 ) n_right++; - if ( strcmp(gas, " 90.00") == 0 ) n_right++; - - /* Hexgonal or monoclinic? */ - if ( n_right == 2 ) { - - if ( (strcmp(als, " 120.00") == 0) - && (strcmp(bs, cs) == 0) ) - { - /* Hexagonal, unique axis a */ - cell_set_lattice_type(cell, L_HEXAGONAL); - cell_set_unique_axis(cell, 'a'); - return; - } - - if ( (strcmp(bes, " 120.00") == 0) - && (strcmp(as, cs) == 0) ) - { - /* Hexagonal, unique axis b */ - cell_set_lattice_type(cell, L_HEXAGONAL); - cell_set_unique_axis(cell, 'b'); - return; - } - - if ( (strcmp(gas, " 120.00") == 0) - && (strcmp(as, bs) == 0) ) - { - /* Hexagonal, unique axis c */ - cell_set_lattice_type(cell, L_HEXAGONAL); - cell_set_unique_axis(cell, 'c'); - return; - } - - if ( strcmp(als, " 90.00") != 0 ) { - /* Monoclinic, unique axis a */ - cell_set_lattice_type(cell, L_MONOCLINIC); - cell_set_unique_axis(cell, 'a'); - return; - } - - if ( strcmp(bes, " 90.00") != 0 ) { - /* Monoclinic, unique axis b */ - cell_set_lattice_type(cell, L_MONOCLINIC); - cell_set_unique_axis(cell, 'b'); - return; - } - - if ( strcmp(gas, " 90.00") != 0 ) { - /* Monoclinic, unique axis c */ - cell_set_lattice_type(cell, L_MONOCLINIC); - cell_set_unique_axis(cell, 'c'); - return; - } - } - - /* Triclinic, unique axis irrelevant. */ - cell_set_lattice_type(cell, L_TRICLINIC); -} - - -UnitCell *load_cell_from_pdb(const char *filename) -{ - FILE *fh; - char *rval; - UnitCell *cell = NULL; - - fh = fopen(filename, "r"); - if ( fh == NULL ) { - ERROR("Couldn't open '%s'\n", filename); - return NULL; - } - - do { - - char line[1024]; - - rval = fgets(line, 1023, fh); - - if ( strncmp(line, "CRYST1", 6) == 0 ) { - - float a, b, c, al, be, ga; - char as[10], bs[10], cs[10]; - char als[8], bes[8], gas[8]; - int r; - - memcpy(as, line+6, 9); as[9] = '\0'; - memcpy(bs, line+15, 9); bs[9] = '\0'; - memcpy(cs, line+24, 9); cs[9] = '\0'; - memcpy(als, line+33, 7); als[7] = '\0'; - memcpy(bes, line+40, 7); bes[7] = '\0'; - memcpy(gas, line+47, 7); gas[7] = '\0'; - - STATUS("'%s' '%s' '%s'\n", as, bs, cs); - STATUS("'%s' '%s' '%s'\n", als, bes, gas); - - r = sscanf(as, "%f", &a); - r += sscanf(bs, "%f", &b); - r += sscanf(cs, "%f", &c); - r += sscanf(als, "%f", &al); - r += sscanf(bes, "%f", &be); - r += sscanf(gas, "%f", &ga); - - if ( r != 6 ) { - STATUS("Couldn't understand CRYST1 line.\n"); - continue; - } - - cell = cell_new_from_parameters(a*1e-10, - b*1e-10, c*1e-10, - deg2rad(al), - deg2rad(be), - deg2rad(ga)); - - determine_lattice(cell, as, bs, cs, als, bes, gas); - - if ( strlen(line) > 65 ) { - cell_set_centering(cell, line[55]); - } else { - cell_set_pointgroup(cell, "1"); - cell_set_spacegroup(cell, "P 1"); - ERROR("CRYST1 line without centering.\n"); - } - - break; /* Done */ - } - - } while ( rval != NULL ); - - fclose(fh); - - validate_cell(cell); - - return cell; -} - - -/* Force the linker to bring in CBLAS to make GSL happy */ -void cell_fudge_gslcblas() -{ - STATUS("%p\n", cblas_sgemm); -} - - -UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot) -{ - UnitCell *out; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double xnew, ynew, znew; - - cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy, - &bsz, &csx, &csy, &csz); - - /* Rotate by "omega" about +z (parallel to c* and c unless triclinic) */ - xnew = asx*cos(omega) + asy*sin(omega); - ynew = -asx*sin(omega) + asy*cos(omega); - znew = asz; - asx = xnew; asy = ynew; asz = znew; - xnew = bsx*cos(omega) + bsy*sin(omega); - ynew = -bsx*sin(omega) + bsy*cos(omega); - znew = bsz; - bsx = xnew; bsy = ynew; bsz = znew; - xnew = csx*cos(omega) + csy*sin(omega); - ynew = -csx*sin(omega) + csy*cos(omega); - znew = csz; - csx = xnew; csy = ynew; csz = znew; - - /* Rotate by "phi" about +x (not parallel to anything specific) */ - xnew = asx; - ynew = asy*cos(phi) + asz*sin(phi); - znew = -asy*sin(phi) + asz*cos(phi); - asx = xnew; asy = ynew; asz = znew; - xnew = bsx; - ynew = bsy*cos(phi) + bsz*sin(phi); - znew = -bsy*sin(phi) + bsz*cos(phi); - bsx = xnew; bsy = ynew; bsz = znew; - xnew = csx; - ynew = csy*cos(phi) + csz*sin(phi); - znew = -csy*sin(phi) + csz*cos(phi); - csx = xnew; csy = ynew; csz = znew; - - /* Rotate by "rot" about the new +z (in-plane rotation) */ - xnew = asx*cos(rot) + asy*sin(rot); - ynew = -asx*sin(rot) + asy*cos(rot); - znew = asz; - asx = xnew; asy = ynew; asz = znew; - xnew = bsx*cos(rot) + bsy*sin(rot); - ynew = -bsx*sin(rot) + bsy*cos(rot); - znew = bsz; - bsx = xnew; bsy = ynew; bsz = znew; - xnew = csx*cos(rot) + csy*sin(rot); - ynew = -csx*sin(rot) + csy*cos(rot); - znew = csz; - csx = xnew; csy = ynew; csz = znew; - - out = cell_new_from_cell(in); - cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); - - return out; -} - - -int cell_is_sensible(UnitCell *cell) -{ - double a, b, c, al, be, ga; - - cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); - if ( al + be + ga >= 2.0*M_PI ) return 0; - if ( al + be - ga >= 2.0*M_PI ) return 0; - if ( al - be + ga >= 2.0*M_PI ) return 0; - if ( - al + be + ga >= 2.0*M_PI ) return 0; - if ( al + be + ga <= 0.0 ) return 0; - if ( al + be - ga <= 0.0 ) return 0; - if ( al - be + ga <= 0.0 ) return 0; - if ( - al + be + ga <= 0.0 ) return 0; - if ( isnan(al) ) return 0; - if ( isnan(be) ) return 0; - if ( isnan(ga) ) return 0; - return 1; -} - - -static int bravais_lattice(UnitCell *cell) -{ - LatticeType lattice = cell_get_lattice_type(cell); - char centering = cell_get_centering(cell); - char ua = cell_get_unique_axis(cell); - - switch ( centering ) - { - case 'P' : - return 1; - - case 'A' : - case 'B' : - case 'C' : - if ( (lattice != L_MONOCLINIC) - && (lattice != L_ORTHORHOMBIC) ) - { - return 0; - } - if ( (ua=='a') && (centering=='A') ) return 1; - if ( (ua=='b') && (centering=='B') ) return 1; - if ( (ua=='c') && (centering=='C') ) return 1; - return 0; - - case 'I' : - if ( (lattice == L_ORTHORHOMBIC) - || (lattice == L_TETRAGONAL) - || (lattice == L_CUBIC) ) - { - return 1; - } - return 0; - - case 'F' : - if ( (lattice == L_ORTHORHOMBIC) || (lattice == L_CUBIC) ) { - return 1; - } - return 0; - - case 'H' : - if ( lattice == L_HEXAGONAL ) return 1; - return 0; - - default : - return 0; - } -} - - -/** - * validate_cell: - * @cell: A %UnitCell to validate - * - * Perform some checks for crystallographic validity @cell, such as that the - * lattice is a conventional Bravais lattice. - * Warnings are printied if any of the checks are failed. - * - */ -void validate_cell(UnitCell *cell) -{ - int err = 0; - - if ( !cell_is_sensible(cell) ) { - ERROR("Warning: Unit cell parameters are not sensible.\n"); - err = 1; - } - - if ( !bravais_lattice(cell) ) { - ERROR("Warning: Unit cell is not a conventional Bravais" - " lattice.\n"); - err = 1; - } - - if ( err ) cell_print(cell); -} diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index 851237cf..e4e9d6ce 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -1,7 +1,7 @@ /* * cell.h * - * Unit Cell Calculations + * A class representing a unit cell * * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. @@ -128,24 +128,6 @@ extern void cell_set_centering(UnitCell *cell, char centering); extern char cell_get_unique_axis(UnitCell *cell); extern void cell_set_unique_axis(UnitCell *cell, char unique_axis); -extern double resolution(UnitCell *cell, - signed int h, signed int k, signed int l); - -extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat); -extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi, - double rot); - -extern void cell_print(UnitCell *cell); - -extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose, - const float *ltl, int reduce); - -extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *tempcell); - -extern UnitCell *load_cell_from_pdb(const char *filename); - -extern int cell_is_sensible(UnitCell *cell); - -extern void validate_cell(UnitCell *cell); +extern const char *cell_rep(UnitCell *cell); #endif /* CELL_H */ diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 24669aef..091b4fed 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -36,6 +36,7 @@ #include "utils.h" #include "cell.h" +#include "cell-utils.h" #include "image.h" #include "peaks.h" #include "beam-parameters.h" diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 102a5312..7d912902 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -49,6 +49,7 @@ #include "index-priv.h" #include "reax.h" #include "geometry.h" +#include "cell-utils.h" /* Base class constructor for unspecialised indexing private data */ diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index b87eb56e..f7f6c650 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -51,6 +51,7 @@ #include "filters.h" #include "reflist-utils.h" #include "beam-parameters.h" +#include "cell-utils.h" /* Degree of polarisation of X-ray beam */ diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c index 5cfa908a..8bd9f29f 100644 --- a/libcrystfel/src/reax.c +++ b/libcrystfel/src/reax.c @@ -47,6 +47,7 @@ #include "utils.h" #include "peaks.h" #include "cell.h" +#include "cell-utils.h" #include "index.h" #include "index-priv.h" diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c index 2e6715f2..f2292929 100644 --- a/libcrystfel/src/reflist-utils.c +++ b/libcrystfel/src/reflist-utils.c @@ -35,6 +35,7 @@ #include "reflist.h" #include "cell.h" +#include "cell-utils.h" #include "utils.h" #include "reflist-utils.h" #include "symmetry.h" -- cgit v1.2.3