diff options
| -rw-r--r-- | libcrystfel/src/cell-utils.c | 340 | ||||
| -rw-r--r-- | libcrystfel/src/cell-utils.h | 6 | ||||
| -rw-r--r-- | libcrystfel/src/reax.c | 4 |
3 files changed, 285 insertions, 65 deletions
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index ca697391..508b70e5 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -39,6 +39,7 @@ #include <gsl/gsl_matrix.h> #include <gsl/gsl_blas.h> #include <gsl/gsl_linalg.h> +#include <assert.h> #include "cell.h" #include "cell-utils.h" @@ -99,6 +100,56 @@ static const char *str_lattice(LatticeType l) } +int right_handed(UnitCell *cell) +{ + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + struct rvec aCb; + double aCb_dot_c; + int rh_reciprocal; + int rh_direct; + + if ( cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz) ) { + ERROR("Couldn't get reciprocal cell.\n"); + return 0; + } + + /* "a" cross "b" */ + aCb.u = asy*bsz - asz*bsy; + aCb.v = - (asx*bsz - asz*bsx); + aCb.w = asx*bsy - asy*bsx; + + /* "a cross b" dot "c" */ + aCb_dot_c = aCb.u*csx + aCb.v*csy + aCb.w*csz; + + rh_reciprocal = aCb_dot_c > 0.0; + + if ( cell_get_cartesian(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz) ) { + ERROR("Couldn't get direct cell.\n"); + return 0; + } + + /* "a" cross "b" */ + aCb.u = asy*bsz - asz*bsy; + aCb.v = - (asx*bsz - asz*bsx); + aCb.w = asx*bsy - asy*bsx; + + /* "a cross b" dot "c" */ + aCb_dot_c = aCb.u*csx + aCb.v*csy + aCb.w*csz; + + rh_direct = aCb_dot_c > 0.0; + + assert(rh_reciprocal == rh_direct); + + return rh_direct; +} + + void cell_print(UnitCell *cell) { double asx, asy, asz; @@ -106,9 +157,31 @@ void cell_print(UnitCell *cell) double csx, csy, csz; double a, b, c, alpha, beta, gamma; double ax, ay, az, bx, by, bz, cx, cy, cz; + LatticeType lt; + + lt = cell_get_lattice_type(cell); - STATUS("%s %c\n", str_lattice(cell_get_lattice_type(cell)), - cell_get_centering(cell)); + STATUS("%s %c", str_lattice(lt), cell_get_centering(cell)); + + switch ( lt ) + { + case L_MONOCLINIC : + case L_TETRAGONAL : + case L_HEXAGONAL : + STATUS(", unique axis %c", cell_get_unique_axis(cell)); + break; + + default : + break; + } + + 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); @@ -134,14 +207,200 @@ void cell_print(UnitCell *cell) 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)); } +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 ) { + if ( (ua=='a') && (centering=='A') ) return 1; + if ( (ua=='b') && (centering=='B') ) return 1; + if ( (ua=='c') && (centering=='C') ) return 1; + } else if ( lattice == L_ORTHORHOMBIC) { + 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; + } +} + + +/* Turn any cell into a primitive one, e.g. for comparison purposes */ +UnitCell *uncenter_cell(UnitCell *in) +{ + UnitCell *cell; + + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + double nax, nay, naz; + double nbx, nby, nbz; + double ncx, ncy, ncz; + const double OT = 1.0/3.0; + const double TT = 2.0/3.0; + const double H = 0.5; + LatticeType lt; + char ua, cen; + + if ( !bravais_lattice(in) ) { + ERROR("Cannot uncenter: not a Bravais lattice.\n"); + return NULL; + } + + cell_get_cartesian(in, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + lt = cell_get_lattice_type(in); + ua = cell_get_unique_axis(in); + cen = cell_get_centering(in); + + if ( ua == 'a' ) { + double tx, ty, tz; + tx = ax; ty = ay; tz = az; + ax = cx; ay = cy; az = cz; + cx = tx; cy = ty; cz = tz; + if ( cen == 'A' ) cen = 'C'; + } + + if ( ua == 'b' ) { + double tx, ty, tz; + tx = bx; ty = by; tz = bz; + ax = cx; ay = cy; az = cz; + cx = tx; cy = ty; cz = tz; + if ( cen == 'B' ) cen = 'C'; + } + + cell = cell_new(); + + switch ( cen ) { + + case 'P' : + case 'R' : + nax = ax; nay = ay; naz = az; + nbx = bx; nby = by; nbz = bz; + ncx = cx; ncy = cy; ncz = cz; /* Identity */ + cell_set_lattice_type(cell, lt); + cell_set_centering(cell, cen); /* P->P, R->R */ + break; + + case 'I' : + nax = - H*ax + H*bx + H*cx; + nay = - H*ay + H*by + H*cy; + naz = - H*az + H*bz + H*cz; + nbx = H*ax - H*bx + H*cx; + nby = H*ay - H*by + H*cy; + nbz = H*az - H*bz + H*cz; + ncx = H*ax + H*bx - H*cx; + ncy = H*ay + H*by - H*cy; + ncz = H*az + H*bz - H*cz; + if ( lt == L_CUBIC ) { + cell_set_lattice_type(cell, L_RHOMBOHEDRAL); + cell_set_centering(cell, 'R'); + } else { + /* Tetragonal or orthorhombic */ + cell_set_lattice_type(cell, L_TRICLINIC); + cell_set_centering(cell, 'P'); + } + break; + + case 'F' : + nax = 0*ax + H*bx + H*cx; + nay = 0*ay + H*by + H*cy; + naz = 0*az + H*bz + H*cz; + nbx = H*ax + 0*bx + H*cx; + nby = H*ay + 0*by + H*cy; + nbz = H*az + 0*bz + H*cz; + ncx = H*ax + H*bx + 0*cx; + ncy = H*ay + H*by + 0*cy; + ncz = H*az + H*bz + 0*cz; + if ( lt == L_CUBIC ) { + cell_set_lattice_type(cell, L_RHOMBOHEDRAL); + cell_set_centering(cell, 'R'); + + } else { + assert(lt == L_ORTHORHOMBIC); + cell_set_lattice_type(cell, L_TRICLINIC); + cell_set_centering(cell, 'P'); + } + break; + + case 'C' : + nax = H*ax + H*bx + 0*cx; + nay = H*ay + H*by + 0*cy; + naz = H*az + H*bz + 0*cz; + nbx = -H*ax + H*bx + 0*cx; + nby = -H*ay + H*by + 0*cy; + nbz = -H*az + H*bz + 0*cz; + ncx = 0*ax + 0*bx + 1*cx; + ncy = 0*ay + 0*by + 1*cy; + ncz = 0*az + 0*bz + 1*cz; + cell_set_lattice_type(cell, L_MONOCLINIC); + cell_set_centering(cell, 'P'); + cell_set_unique_axis(cell, 'c'); + break; + + case 'H' : + nax = TT*ax + OT*bx + OT*cx; + nay = TT*ay + OT*by + OT*cy; + naz = TT*az + OT*bz + OT*cz; + nbx = - OT*ax + OT*bx + OT*cx; + nby = - OT*ay + OT*by + OT*cy; + nbz = - OT*az + OT*bz + OT*cz; + ncx = - OT*ax - TT*bx + OT*cx; + ncy = - OT*ay - TT*by + OT*cy; + ncz = - OT*az - TT*bz + OT*cz; + assert(lt == L_HEXAGONAL); + cell_set_lattice_type(cell, L_RHOMBOHEDRAL); + cell_set_centering(cell, 'R'); + break; + + default : + ERROR("Invalid centering '%c'\n", cell_get_centering(in)); + return NULL; + + } + + cell_set_cartesian(cell, nax, nay, naz, nbx, nby, nbz, ncx, ncy, ncz); + + return cell; +} + + #define MAX_CAND (1024) -static int right_handed(struct rvec a, struct rvec b, struct rvec c) +static int right_handed_vec(struct rvec a, struct rvec b, struct rvec c) { struct rvec aCb; double aCb_dot_c; @@ -178,7 +437,7 @@ static int same_vector(struct cvec a, struct cvec b) /* Attempt to make 'cell' fit into 'template' somehow */ -UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, +UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, const float *tols, int reduce) { signed int n1l, n2l, n3l; @@ -194,6 +453,11 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, int ncand[3] = {0,0,0}; signed int ilow, ihigh; float angtol = deg2rad(tols[3]); + UnitCell *cell; + UnitCell *template; + + cell = uncenter_cell(cell_in); + template = uncenter_cell(template_in); if ( cell_get_reciprocal(template, &asx, &asy, &asz, &bsx, &bsy, &bsz, @@ -350,8 +614,8 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, 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; + if ( !right_handed_vec(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 @@ -464,7 +728,7 @@ UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template) } /* Flip c if not right-handed */ - if ( !right_handed(real_a, real_b, real_c) ) { + if ( !right_handed_vec(real_a, real_b, real_c) ) { real_c.u = -real_c.u; real_c.v = -real_c.v; real_c.w = -real_c.w; @@ -690,12 +954,6 @@ UnitCell *load_cell_from_pdb(const char *filename) fclose(fh); - /* FIXME: Turn "H" centered cells into "R" cells */ - if ( cell_get_centering(cell) == 'H' ) { - STATUS("Turning your H-centered (PDB convention) cell into" - " an R-centered one.\n"); - } - validate_cell(cell); return cell; @@ -789,55 +1047,6 @@ int cell_is_sensible(UnitCell *cell) } -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 @@ -862,5 +1071,10 @@ void validate_cell(UnitCell *cell) err = 1; } + if ( !right_handed(cell) ) { + ERROR("Warning: Unit cell is not right handed.\n"); + err = 1; + } + if ( err ) cell_print(cell); } diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index 8d8cfc6d..b2cb7b67 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -55,4 +55,10 @@ extern int cell_is_sensible(UnitCell *cell); extern void validate_cell(UnitCell *cell); +extern UnitCell *uncenter_cell(UnitCell *in); + +extern int bravais_lattice(UnitCell *cell); + +extern int right_handed(UnitCell *cell); + #endif /* CELL_UTILS_H */ diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c index 8bd9f29f..29ea7294 100644 --- a/libcrystfel/src/reax.c +++ b/libcrystfel/src/reax.c @@ -745,7 +745,7 @@ static double max_feature_resolution(ImageFeatureList *flist) } -static int right_handed(struct rvec a, struct rvec b, struct rvec c) +static int right_handed_vec(struct rvec a, struct rvec b, struct rvec c) { struct rvec aCb; double aCb_dot_c; @@ -958,7 +958,7 @@ static void assemble_cells_from_candidates(struct image *image, bi.u = vj.x; bi.v = vj.y; bi.w = vj.z; ci.u = vk.x; ci.v = vk.y; ci.w = vk.z; - if ( !right_handed(ai, bi, ci) ) continue; + if ( !right_handed_vec(ai, bi, ci) ) continue; /* We have three vectors with the right angles */ cnew = cell_new_from_direct_axes(ai, bi, ci); |