From 766d90e7d92ee0e10877e90466829ff3f4925269 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 1 Jun 2012 18:15:17 +0200 Subject: Handle lattice type, centering and unique axis information --- libcrystfel/src/cell.c | 280 +++++++++++++++++++++++++++++++++++++++++++++---- libcrystfel/src/cell.h | 22 ++++ 2 files changed, 281 insertions(+), 21 deletions(-) diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 3d2a28bf..1cd6fb88 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -92,8 +92,11 @@ struct _unitcell { double ays; double bys; double cys; double azs; double bzs; double czs; - char *pointgroup; - char *spacegroup; + char *pointgroup; + char *spacegroup; + LatticeType lattice_type; + char centering; + char unique_axis; }; @@ -126,6 +129,9 @@ UnitCell *cell_new() cell->pointgroup = strdup("1"); cell->spacegroup = strdup("P 1"); + cell->lattice_type = L_TRICLINIC; + cell->centering = 'P'; + cell->unique_axis = 'c'; return cell; } @@ -296,6 +302,24 @@ void cell_set_pointgroup(UnitCell *cell, const char *sym) } +void cell_set_centering(UnitCell *cell, char centering) +{ + cell->centering = centering; +} + + +void cell_set_lattice_type(UnitCell *cell, LatticeType lattice_type) +{ + cell->lattice_type = lattice_type; +} + + +void cell_set_unique_axis(UnitCell *cell, char unique_axis) +{ + cell->unique_axis = unique_axis; +} + + /************************* Getter helper functions ****************************/ static int cell_crystallographic_to_cartesian(UnitCell *cell, @@ -561,6 +585,23 @@ const char *cell_get_spacegroup(UnitCell *cell) } +char cell_get_centering(UnitCell *cell) +{ + return cell->centering; +} + + +LatticeType cell_get_lattice_type(UnitCell *cell) +{ + return cell->lattice_type; +} + + +char cell_get_unique_axis(UnitCell *cell) +{ + return cell->unique_axis; +} + @@ -617,6 +658,23 @@ UnitCell *cell_rotate(UnitCell *in, struct quaternion quat) } +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; @@ -625,6 +683,9 @@ void cell_print(UnitCell *cell) 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"); @@ -1019,21 +1080,122 @@ double resolution(UnitCell *cell, signed int h, signed int k, signed int l) } -static void cell_set_pointgroup_from_pdb(UnitCell *cell, const char *sym) +static void determine_lattice(UnitCell *cell, + const char *as, const char *bs, const char *cs, + const char *als, const char *bes, const char *gas) { - char *new = NULL; + int n_right; - if ( strcmp(sym, "P 1") == 0 ) new = "1"; - if ( strcmp(sym, "P 63") == 0 ) new = "6"; - if ( strcmp(sym, "P 21 21 21") == 0 ) new = "222"; - if ( strcmp(sym, "P 2 2 2") == 0 ) new = "222"; - if ( strcmp(sym, "P 43 21 2") == 0 ) new = "422"; + /* Rhombohedral or cubic? */ + if ( (strcmp(as, bs) == 0) && (strcmp(as, cs) == 0) ) { - if ( new != NULL ) { - cell_set_pointgroup(cell, new); - } else { - ERROR("Can't determine point group for '%s'\n", sym); + 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); } @@ -1060,7 +1222,6 @@ UnitCell *load_cell_from_pdb(const char *filename) float a, b, c, al, be, ga; char as[10], bs[10], cs[10]; char als[8], bes[8], gas[8]; - char *sym; int r; memcpy(as, line+6, 9); as[9] = '\0'; @@ -1088,16 +1249,14 @@ UnitCell *load_cell_from_pdb(const char *filename) deg2rad(be), deg2rad(ga)); + determine_lattice(cell, as, bs, cs, als, bes, gas); + if ( strlen(line) > 65 ) { - sym = strndup(line+55, 10); - notrail(sym); - cell_set_pointgroup_from_pdb(cell, sym); - cell_set_spacegroup(cell, sym); - free(sym); + cell_set_centering(cell, line[55]); } else { - cell_set_pointgroup_from_pdb(cell, "P 1"); + cell_set_pointgroup(cell, "1"); cell_set_spacegroup(cell, "P 1"); - ERROR("CRYST1 line without space group.\n"); + ERROR("CRYST1 line without centering.\n"); } break; /* Done */ @@ -1107,6 +1266,8 @@ UnitCell *load_cell_from_pdb(const char *filename) fclose(fh); + validate_cell(cell); + return cell; } @@ -1196,3 +1357,80 @@ int cell_is_sensible(UnitCell *cell) 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 bd2719dd..851237cf 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -48,6 +48,17 @@ struct rvec double w; }; +typedef enum +{ + L_TRICLINIC, + L_MONOCLINIC, + L_ORTHORHOMBIC, + L_TETRAGONAL, + L_RHOMBOHEDRAL, + L_HEXAGONAL, + L_CUBIC +} LatticeType; + /** * UnitCell: @@ -108,6 +119,15 @@ extern const char *cell_get_pointgroup(UnitCell *cell); extern const char *cell_get_spacegroup(UnitCell *cell); +extern LatticeType cell_get_lattice_type(UnitCell *cell); +extern void cell_set_lattice_type(UnitCell *cell, LatticeType lattice_type); + +extern char cell_get_centering(UnitCell *cell); +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); @@ -126,4 +146,6 @@ extern UnitCell *load_cell_from_pdb(const char *filename); extern int cell_is_sensible(UnitCell *cell); +extern void validate_cell(UnitCell *cell); + #endif /* CELL_H */ -- cgit v1.2.3 From 32fdb97adf4f03e0c3ddceb96e4fe6111c4125c0 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Jun 2012 10:57:11 +0200 Subject: WIP --- libcrystfel/src/cell.c | 3 +++ libcrystfel/src/mosflm.c | 57 ++++++++++++++++++++++++++++++++++++++---------- 2 files changed, 48 insertions(+), 12 deletions(-) diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 1cd6fb88..8f380128 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -1231,6 +1231,9 @@ UnitCell *load_cell_from_pdb(const char *filename) 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); diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index ed118aa4..a2d8520e 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -272,14 +272,54 @@ static void mosflm_sendline(const char *line, struct mosflm_data *mosflm) } +static const char *spacegroup_for_lattice(UnitCell *cell) +{ + LatticeType latt; + char centering; + + latt = cell_get_lattice_type(cell); + centering = cell_get_centering(cell); + + switch ( latt ) + { + case L_TRICLINIC : + g = "1"; + break; + + case L_MONOCLINIC : + g = "2"; + break; + + case L_ORTHORHOMBIC : + g = "222"; + break; + + case L_TETRAGONAL : + g = "4"; + break; + + case L_RHOMBOHEDRAL : + g = "3"; + break; + + case L_HEXAGONAL : + g = "6"; + break; + + case L_CUBIC : + g = "23"; + break; + + + return "P1"; +} + + static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) { char tmp[256]; - char symm[32]; - const char *sg; double wavelength; double a, b, c, alpha, beta, gamma; - int i, j; switch ( mosflm->step ) { @@ -305,15 +345,8 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) case 3 : if ( mosflm->target_cell != NULL ) { - sg = cell_get_spacegroup(mosflm->target_cell); - /* Remove white space from space group */ - j = 0; - for ( i=0; itarget_cell); snprintf(tmp, 255, "SYMM %s\n", symm); mosflm_sendline(tmp, mosflm); } else { -- cgit v1.2.3 From 38f4df68678bb4c6eb8cbb3b721db2e130643299 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 7 Aug 2012 17:38:33 +0200 Subject: "divisors" not needed any more Was it ever needed? --- libcrystfel/src/symmetry.c | 3 --- 1 file changed, 3 deletions(-) diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 0352e441..68b39bcb 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -68,7 +68,6 @@ struct _symoplist int n_ops; int max_ops; char *name; - int *divisors; int num_equivs; }; @@ -84,7 +83,6 @@ struct _symopmask static void alloc_ops(SymOpList *ops) { ops->ops = realloc(ops->ops, ops->max_ops*sizeof(struct sym_op)); - ops->divisors = realloc(ops->divisors, ops->max_ops*sizeof(int)); } @@ -127,7 +125,6 @@ static SymOpList *new_symoplist() new->max_ops = 16; new->n_ops = 0; new->ops = NULL; - new->divisors = NULL; new->name = NULL; new->num_equivs = 1; alloc_ops(new); -- cgit v1.2.3 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 + src/check_hkl.c | 1 + src/compare_hkl.c | 1 + src/indexamajig.c | 1 + src/partial_sim.c | 17 +- src/pattern_sim.c | 1 + src/post-refinement.c | 1 + src/powder_plot.c | 1 + src/render_hkl.c | 1 + tests/gpu_sim_check.c | 1 + tests/pr_gradient_check.c | 1 + 20 files changed, 951 insertions(+), 853 deletions(-) create mode 100644 libcrystfel/src/cell-utils.c create mode 100644 libcrystfel/src/cell-utils.h 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" diff --git a/src/check_hkl.c b/src/check_hkl.c index 1dfe91a0..84cfdccf 100644 --- a/src/check_hkl.c +++ b/src/check_hkl.c @@ -43,6 +43,7 @@ #include "symmetry.h" #include "reflist.h" #include "reflist-utils.h" +#include "cell-utils.h" /* Number of bins for plot of resolution shells */ diff --git a/src/compare_hkl.c b/src/compare_hkl.c index b6a525fc..2b0e5f90 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -44,6 +44,7 @@ #include "statistics.h" #include "symmetry.h" #include "reflist-utils.h" +#include "cell-utils.h" /* Number of bins for plot of resolution shells */ diff --git a/src/indexamajig.c b/src/indexamajig.c index 074be20a..4b3dcb5d 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -62,6 +62,7 @@ #include "geometry.h" #include "stream.h" #include "reflist-utils.h" +#include "cell-utils.h" #include "im-sandbox.h" diff --git a/src/partial_sim.c b/src/partial_sim.c index 62c72ceb..e9a7003e 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -40,14 +40,15 @@ #include #include -#include -#include -#include -#include -#include -#include -#include -#include +#include "image.h" +#include "utils.h" +#include "reflist-utils.h" +#include "symmetry.h" +#include "beam-parameters.h" +#include "geometry.h" +#include "stream.h" +#include "thread-pool.h" +#include "cell-utils.h" /* Number of bins for partiality graph */ #define NBINS 50 diff --git a/src/pattern_sim.c b/src/pattern_sim.c index f9831e79..06a1f5e5 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -42,6 +42,7 @@ #include "diffraction.h" #include "diffraction-gpu.h" #include "cell.h" +#include "cell-utils.h" #include "utils.h" #include "hdf5-file.h" #include "detector.h" diff --git a/src/post-refinement.c b/src/post-refinement.c index 78c47d07..ac4de1bc 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -45,6 +45,7 @@ #include "symmetry.h" #include "geometry.h" #include "cell.h" +#include "cell-utils.h" /* Maximum number of iterations of NLSq to do for each image per macrocycle. */ diff --git a/src/powder_plot.c b/src/powder_plot.c index dd7f3f29..7917ccd0 100644 --- a/src/powder_plot.c +++ b/src/powder_plot.c @@ -51,6 +51,7 @@ #include "beam-parameters.h" #include "reflist-utils.h" #include "symmetry.h" +#include "cell-utils.h" struct bin_stats { diff --git a/src/render_hkl.c b/src/render_hkl.c index 9e7a6dda..3ae3b345 100644 --- a/src/render_hkl.c +++ b/src/render_hkl.c @@ -48,6 +48,7 @@ #include "render_hkl.h" #include "reflist.h" #include "reflist-utils.h" +#include "cell-utils.h" #define KEY_FILENAME "key.pdf" diff --git a/tests/gpu_sim_check.c b/tests/gpu_sim_check.c index 5b6f6da7..5c45ea02 100644 --- a/tests/gpu_sim_check.c +++ b/tests/gpu_sim_check.c @@ -36,6 +36,7 @@ #include #include #include +#include #ifdef HAVE_CLOCK_GETTIME diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index 6972b66a..b39dfd05 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -32,6 +32,7 @@ #include #include +#include #include #include #include "../src/post-refinement.h" -- cgit v1.2.3 From 49e99dd17d003df568ef676a2c5b66ed82193dc7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 7 Aug 2012 18:20:57 +0200 Subject: Remove space group from unit cell --- libcrystfel/src/cell.c | 17 ----------------- libcrystfel/src/cell.h | 3 --- 2 files changed, 20 deletions(-) diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 64bcd3f4..37e201e3 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -89,7 +89,6 @@ struct _unitcell { double azs; double bzs; double czs; char *pointgroup; - char *spacegroup; LatticeType lattice_type; char centering; char unique_axis; @@ -124,7 +123,6 @@ UnitCell *cell_new() cell->rep = CELL_REP_CRYST; cell->pointgroup = strdup("1"); - cell->spacegroup = strdup("P 1"); cell->lattice_type = L_TRICLINIC; cell->centering = 'P'; cell->unique_axis = 'c'; @@ -144,7 +142,6 @@ void cell_free(UnitCell *cell) { if ( cell == NULL ) return; free(cell->pointgroup); - free(cell->spacegroup); free(cell); } @@ -263,7 +260,6 @@ UnitCell *cell_new_from_cell(UnitCell *orig) cell_get_cartesian(orig, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); cell_set_cartesian(new, ax, ay, az, bx, by, bz, cx, cy, cz); cell_set_pointgroup(new, orig->pointgroup); - cell_set_spacegroup(new, orig->spacegroup); return new; } @@ -284,13 +280,6 @@ void cell_set_reciprocal(UnitCell *cell, } -void cell_set_spacegroup(UnitCell *cell, const char *sym) -{ - free(cell->spacegroup); - cell->spacegroup = strdup(sym); -} - - void cell_set_pointgroup(UnitCell *cell, const char *sym) { free(cell->pointgroup); @@ -575,12 +564,6 @@ const char *cell_get_pointgroup(UnitCell *cell) } -const char *cell_get_spacegroup(UnitCell *cell) -{ - return cell->spacegroup; -} - - char cell_get_centering(UnitCell *cell) { return cell->centering; diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index e4e9d6ce..4e90a8ad 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -93,7 +93,6 @@ extern void cell_set_parameters(UnitCell *cell, double a, double b, double c, extern void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az); extern void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz); extern void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz); -extern void cell_set_spacegroup(UnitCell *cell, const char *sym); extern void cell_set_pointgroup(UnitCell *cell, const char *sym); @@ -117,8 +116,6 @@ extern void cell_set_reciprocal(UnitCell *cell, extern const char *cell_get_pointgroup(UnitCell *cell); -extern const char *cell_get_spacegroup(UnitCell *cell); - extern LatticeType cell_get_lattice_type(UnitCell *cell); extern void cell_set_lattice_type(UnitCell *cell, LatticeType lattice_type); -- cgit v1.2.3 From e67a04b7492806816f2c5ecbadd9814eb33068e1 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 8 Aug 2012 11:20:14 +0200 Subject: WIP --- libcrystfel/src/cell-utils.c | 5 ++++- libcrystfel/src/mosflm.c | 14 +++++++++++--- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index c30d9a03..ca697391 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -378,7 +378,6 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, } - UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template) { double ax, ay, az; @@ -692,6 +691,10 @@ 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); diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index a2d8520e..e81f20b6 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -276,6 +276,8 @@ static const char *spacegroup_for_lattice(UnitCell *cell) { LatticeType latt; char centering; + char *g = NULL; + char *result; latt = cell_get_lattice_type(cell); centering = cell_get_centering(cell); @@ -291,7 +293,7 @@ static const char *spacegroup_for_lattice(UnitCell *cell) break; case L_ORTHORHOMBIC : - g = "222"; + g = "2 2 2"; break; case L_TETRAGONAL : @@ -307,11 +309,17 @@ static const char *spacegroup_for_lattice(UnitCell *cell) break; case L_CUBIC : - g = "23"; + g = "2 3"; break; + } + assert(g != NULL); + + result = malloc(32); + if ( result == NULL ) return NULL; + snprintf(result, 31, "%c %s", centering, g); - return "P1"; + return result; } -- cgit v1.2.3 From 7a7202862defd2b3638fbe16b7ec315c6c871321 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 29 Aug 2012 16:30:36 +0200 Subject: Add tests/centering_check --- Makefile.am | 9 ++-- tests/.gitignore | 1 + tests/centering_check.c | 113 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 120 insertions(+), 3 deletions(-) create mode 100644 tests/centering_check.c diff --git a/Makefile.am b/Makefile.am index 44e707ea..953a04df 100644 --- a/Makefile.am +++ b/Makefile.am @@ -7,12 +7,13 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \ src/partialator src/check_hkl src/partial_sim noinst_PROGRAMS = tests/list_check tests/integration_check \ - tests/pr_gradient_check tests/symmetry_check + tests/pr_gradient_check tests/symmetry_check \ + tests/centering_check TESTS = tests/list_check tests/first_merge_check tests/second_merge_check \ tests/third_merge_check tests/fourth_merge_check \ tests/integration_check tests/pr_gradient_check \ - tests/symmetry_check + tests/symmetry_check tests/centering_check EXTRA_DIST += tests/first_merge_check tests/second_merge_check \ tests/third_merge_check tests/fourth_merge_check @@ -80,10 +81,12 @@ tests_integration_check_SOURCES = tests/integration_check.c tests_symmetry_check_SOURCES = tests/symmetry_check.c - tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c \ src/post-refinement.c +tests_centering_check_SOURCES = tests/centering_check.c libcrystfel/src/cell.c \ + libcrystfel/src/cell-utils.c + INCLUDES = -I$(top_srcdir)/libcrystfel/src -I$(top_srcdir)/data EXTRA_DIST += src/dw-hdfsee.h src/hdfsee.h src/render_hkl.h \ diff --git a/tests/.gitignore b/tests/.gitignore index c2a45f5a..16acffc3 100644 --- a/tests/.gitignore +++ b/tests/.gitignore @@ -5,4 +5,5 @@ gpu_sim_check integration_check pr_gradient_check symmetry_check +centering_check .dirstamp diff --git a/tests/centering_check.c b/tests/centering_check.c new file mode 100644 index 00000000..0dc79a5f --- /dev/null +++ b/tests/centering_check.c @@ -0,0 +1,113 @@ +/* + * centering_check.c + * + * Check that centering of cells works + * + * Copyright © 2012 Thomas White + * + * 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 + + +static int check_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 ( !right_handed(cell) ) { + ERROR("Warning: Unit cell is not right handed.\n"); + err = 1; + } + + return err; +} + + +static int check_centering() +{ + UnitCell *cell; + UnitCell *n; + int fail = 0; + + cell = cell_new(); + + cell_set_parameters(cell, 10e-10, 20e-10, 30e-10, + deg2rad(90.0), deg2rad(90.0), deg2rad(100.0)); + cell_set_centering(cell, 'C'); + cell_set_lattice_type(cell, L_MONOCLINIC); + cell_set_unique_axis(cell, 'c'); + if ( check_cell(cell) ) fail = 1; + n = uncenter_cell(cell); + if ( check_cell(n) ) fail = 1; + + return fail; +} + + + +int main(int argc, char *argv[]) +{ + int fail = 0; + + /* Triclinic P */ + check_centering(10e-10, 20e-10, 30e-10, 100.0, 120.0, 140.0, + L_TRICLINIC, 'P', '*', + ); + /* Monoclinic P */ + /* Monoclinic A */ + /* Monoclinic B */ + /* Monoclinic C */ + /* Orthorhombic P */ + /* Orthorhombic A */ + /* Orthorhombic B */ + /* Orthorhombic C */ + /* Orthorhombic I */ + /* Orthorhombic F */ + /* Tetragonal P */ + /* Tetragonal I */ + /* Rhombohedral R */ + /* Hexagonal P */ + /* Hexagonal H (PDB-speak for rhombohedral) */ + /* Cubic P */ + /* Cubic I */ + /* Cubic F */ + + + return fail; +} -- cgit v1.2.3 From cf1ed4e4fb4004e89f2e307931a1de4c024e13de Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 29 Aug 2012 16:31:59 +0200 Subject: Send setting information to MOSFLM --- libcrystfel/src/mosflm.c | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index e81f20b6..d5356bab 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -272,15 +272,19 @@ static void mosflm_sendline(const char *line, struct mosflm_data *mosflm) } +/* Turn what we know about the unit cell into something which we can give to + * MOSFLM to make it give us only indexing results compatible with the cell. */ static const char *spacegroup_for_lattice(UnitCell *cell) { LatticeType latt; char centering; + char ua; char *g = NULL; char *result; latt = cell_get_lattice_type(cell); centering = cell_get_centering(cell); + ua = cell_get_unique_axis(cell); switch ( latt ) { @@ -289,7 +293,10 @@ static const char *spacegroup_for_lattice(UnitCell *cell) break; case L_MONOCLINIC : - g = "2"; + /* "2 1 1", "1 2 1" or "1 1 2" depending on unique axis */ + if ( ua == 'a' ) g = "2 1 1"; + if ( ua == 'b' ) g = "1 2 1"; + if ( ua == 'c' ) g = "1 1 2"; break; case L_ORTHORHOMBIC : @@ -297,7 +304,10 @@ static const char *spacegroup_for_lattice(UnitCell *cell) break; case L_TETRAGONAL : - g = "4"; + /* "4 1 1", "1 4 1" or "1 1 4" depending on unique axis */ + if ( ua == 'a' ) g = "4 1 1"; + if ( ua == 'b' ) g = "1 4 1"; + if ( ua == 'c' ) g = "1 1 4"; break; case L_RHOMBOHEDRAL : @@ -305,7 +315,10 @@ static const char *spacegroup_for_lattice(UnitCell *cell) break; case L_HEXAGONAL : - g = "6"; + /* "6 1 1", "1 6 1" or "1 1 6" depending on unique axis */ + if ( ua == 'a' ) g = "6 1 1"; + if ( ua == 'b' ) g = "1 6 1"; + if ( ua == 'c' ) g = "6"; break; case L_CUBIC : @@ -317,7 +330,7 @@ static const char *spacegroup_for_lattice(UnitCell *cell) result = malloc(32); if ( result == NULL ) return NULL; - snprintf(result, 31, "%c %s", centering, g); + snprintf(result, 31, "%c%s", centering, g); return result; } -- cgit v1.2.3 From 2f5af19403e4e13e36d61e85a9dcf77d8c0d82cd Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 29 Aug 2012 16:32:20 +0200 Subject: Don't send cell parameters to MOSFLM --- libcrystfel/src/mosflm.c | 29 +++++++---------------------- 1 file changed, 7 insertions(+), 22 deletions(-) diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index d5356bab..63166919 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -340,7 +340,6 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) { char tmp[256]; double wavelength; - double a, b, c, alpha, beta, gamma; switch ( mosflm->step ) { @@ -351,20 +350,6 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) break; case 2 : - if ( mosflm->target_cell != NULL ) { - cell_get_parameters(mosflm->target_cell, &a, &b, &c, - &alpha, &beta, &gamma); - snprintf(tmp, 255, - "CELL %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n", - a*1e10, b*1e10, c*1e10, - rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); - mosflm_sendline(tmp, mosflm); - } else { - mosflm_sendline("# Do nothing\n", mosflm); - } - break; - - case 3 : if ( mosflm->target_cell != NULL ) { const char *symm; symm = spacegroup_for_lattice(mosflm->target_cell); @@ -375,31 +360,31 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) } break; - case 4 : + case 3 : mosflm_sendline("DISTANCE 67.8\n", mosflm); break; - case 5 : + case 4 : mosflm_sendline("BEAM 0.0 0.0\n", mosflm); break; - case 6 : + case 5 : wavelength = image->lambda*1e10; snprintf(tmp, 255, "WAVELENGTH %10.5f\n", wavelength); mosflm_sendline(tmp, mosflm); break; - case 7 : + case 6 : snprintf(tmp, 255, "NEWMAT %s\n", mosflm->newmatfile); mosflm_sendline(tmp, mosflm); break; - case 8 : + case 7 : snprintf(tmp, 255, "IMAGE %s phi 0 0\n", mosflm->imagefile); mosflm_sendline(tmp, mosflm); break; - case 9 : + case 8 : snprintf(tmp, 255, "AUTOINDEX DPS FILE %s" " IMAGE 1 MAXCELL 1000 REFINE\n", mosflm->sptfile); @@ -411,7 +396,7 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) mosflm_sendline(tmp, mosflm); break; - case 10 : + case 9 : mosflm_sendline("GO\n", mosflm); mosflm->finished_ok = 1; break; -- cgit v1.2.3 From 60ec4009e4bc28ab9ed772ee6fcd8c80c533dccd Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 29 Aug 2012 16:33:21 +0200 Subject: Uncenter the cell before using it for indexing stuff --- libcrystfel/src/cell-utils.c | 340 +++++++++++++++++++++++++++++++++++-------- libcrystfel/src/cell-utils.h | 6 + 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 #include #include +#include #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); -- cgit v1.2.3 From 4b3ffa51ec169406185a76016a29833bc9637264 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 31 Aug 2012 16:20:21 +0200 Subject: WIP on cell transformations --- libcrystfel/src/cell-utils.c | 81 +++++++++++++++++++-------- libcrystfel/src/cell-utils.h | 4 +- libcrystfel/src/cell.c | 54 ++++++++++++++++++ libcrystfel/src/cell.h | 15 +++++ tests/centering_check.c | 127 ++++++++++++++++++++++++++++++++++++------- 5 files changed, 239 insertions(+), 42 deletions(-) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 508b70e5..776526aa 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -83,7 +83,7 @@ UnitCell *cell_rotate(UnitCell *in, struct quaternion quat) } -static const char *str_lattice(LatticeType l) +const char *str_lattice(LatticeType l) { switch ( l ) { @@ -253,17 +253,19 @@ int bravais_lattice(UnitCell *cell) if ( lattice == L_HEXAGONAL ) return 1; return 0; + case 'R' : + if ( lattice == L_RHOMBOHEDRAL ) return 1; + return 0; + default : return 0; } } -/* Turn any cell into a primitive one, e.g. for comparison purposes */ -UnitCell *uncenter_cell(UnitCell *in) +static UnitCellTransformation *uncentering_transformation(UnitCell *in) { - UnitCell *cell; - + UnitCellTransformation *t; double ax, ay, az; double bx, by, bz; double cx, cy, cz; @@ -276,34 +278,29 @@ UnitCell *uncenter_cell(UnitCell *in) 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' ) { + if ( (ua == 'a') || (cen == 'A') ) { double tx, ty, tz; tx = ax; ty = ay; tz = az; ax = cx; ay = cy; az = cz; - cx = tx; cy = ty; cz = tz; + cx = -tx; cy = -ty; cz = -tz; if ( cen == 'A' ) cen = 'C'; + ua = 'c'; } - if ( ua == 'b' ) { + if ( (ua == 'b') || (cen == 'B') ) { double tx, ty, tz; tx = bx; ty = by; tz = bz; - ax = cx; ay = cy; az = cz; - cx = tx; cy = ty; cz = tz; + bx = cx; by = cy; bz = cz; + cx = -tx; cy = -ty; cz = -tz; if ( cen == 'B' ) cen = 'C'; + ua = 'c'; } - cell = cell_new(); - switch ( cen ) { case 'P' : @@ -392,9 +389,39 @@ UnitCell *uncenter_cell(UnitCell *in) } - cell_set_cartesian(cell, nax, nay, naz, nbx, nby, nbz, ncx, ncy, ncz); - return cell; +} + + +/** + * uncenter_cell: + * @in: A %UnitCell + * @t: Location at which to store the transformation which was used. + * + * Turns any cell into a primitive one, e.g. for comparison purposes. The + * transformation which was used is stored at @t, which can be NULL if the + * transformation is not required. + * + * Returns: a primitive version of @in in a conventional (unique axis c) + * setting. + * + */ +UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t) +{ + UnitCell *cell; + UnitCellTransformation *tt; + + if ( !bravais_lattice(in) ) { + ERROR("Cannot uncenter: not a Bravais lattice.\n"); + return NULL; + } + + tt = uncentering_transformation(in); + if ( tt == NULL ) return NULL; + + if ( t != NULL ) *t = tt; + + return cell_transform(in, tt); } @@ -455,9 +482,15 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, float angtol = deg2rad(tols[3]); UnitCell *cell; UnitCell *template; + UnitCellTransformation *to_given_cell; + UnitCell *new_cell_trans; - cell = uncenter_cell(cell_in); - template = uncenter_cell(template_in); + /* "Un-center" the template unit cell to make the comparison easier */ + template = uncenter_cell(template_in, &to_given_cell); + + /* The candidate cell is also uncentered, because it might be centered + * if it came from (e.g.) MOSFLM */ + cell = uncenter_cell(cell_in, NULL); if ( cell_get_reciprocal(template, &asx, &asy, &asz, &bsx, &bsy, &bsz, @@ -638,7 +671,11 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, free(cand[1]); free(cand[2]); - return new_cell; + /* Reverse the de-centering transformation */ + new_cell_trans = cell_transform_inverse(new_cell, to_given_cell); + cell_free(new_cell); + + return new_cell_trans; } diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index b2cb7b67..8acf2f85 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -55,10 +55,12 @@ extern int cell_is_sensible(UnitCell *cell); extern void validate_cell(UnitCell *cell); -extern UnitCell *uncenter_cell(UnitCell *in); +extern UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **tr); extern int bravais_lattice(UnitCell *cell); extern int right_handed(UnitCell *cell); +extern const char *str_lattice(LatticeType l); + #endif /* CELL_UTILS_H */ diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 37e201e3..e2617dfc 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -599,3 +599,57 @@ const char *cell_rep(UnitCell *cell) return "unknown"; } + + +static UnitCellTransformation *inverse_transformation(UnitCellTransformation *t) +{ + /* FIXME: Implementation */ + return NULL; +} + + +/** + * cell_transform: + * @cell: A %UnitCell. + * @t: A %UnitCellTransformation. + * + * Applies @t to @cell. + * + * Returns: Transformed copy of @cell. + * + */ +UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t) +{ + UnitCell *out; + + if ( t == NULL ) return NULL; + + out = cell_new(); + if ( out == NULL ) return NULL; + + /* FIXME: Implementation */ + + return out; +} + + +/** + * cell_transform: + * @cell: A %UnitCell. + * @t: A %UnitCellTransformation. + * + * Applies the inverse of @t to @cell. + * + * Returns: Transformed copy of @cell. + * + */ +UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t) +{ + return cell_transform(cell, inverse_transformation(t)); +} + + +void cell_transformation_print(UnitCellTransformation *t) +{ + +} diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index 4e90a8ad..625884e1 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -68,6 +68,15 @@ typedef enum **/ typedef struct _unitcell UnitCell; + +/** + * UnitCellTransformation: + * + * This opaque data structure represents a tranformation of a unit cell, such + * as a rotation or a centering operation. + **/ +typedef struct _unitcelltransformation UnitCellTransformation; + extern UnitCell *cell_new(void); extern UnitCell *cell_new_from_cell(UnitCell *orig); extern void cell_free(UnitCell *cell); @@ -127,4 +136,10 @@ extern void cell_set_unique_axis(UnitCell *cell, char unique_axis); extern const char *cell_rep(UnitCell *cell); +extern UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t); +extern UnitCell *cell_transform_inverse(UnitCell *cell, + UnitCellTransformation *t); + +extern void cell_transformation_print(UnitCellTransformation *t); + #endif /* CELL_H */ diff --git a/tests/centering_check.c b/tests/centering_check.c index 0dc79a5f..96326d21 100644 --- a/tests/centering_check.c +++ b/tests/centering_check.c @@ -35,46 +35,63 @@ #include -static int check_cell(UnitCell *cell) +static int check_cell(UnitCell *cell, const char *text) { int err = 0; if ( !cell_is_sensible(cell) ) { - ERROR("Warning: Unit cell parameters are not sensible.\n"); + ERROR(" %s unit cell parameters are not sensible.\n", text); err = 1; } if ( !bravais_lattice(cell) ) { - ERROR("Warning: Unit cell is not a conventional Bravais" - " lattice.\n"); + ERROR(" %s unit cell is not a conventional Bravais" + " lattice.\n", text); err = 1; } if ( !right_handed(cell) ) { - ERROR("Warning: Unit cell is not right handed.\n"); + ERROR(" %s unit cell is not right handed.\n", text); err = 1; } + if ( err ) cell_print(cell); + return err; } -static int check_centering() +static int check_centering(double a, double b, double c, + double al, double be, double ga, + LatticeType latt, char cen, char ua) { UnitCell *cell; UnitCell *n; + UnitCellTransformation *t; int fail = 0; - cell = cell_new(); + STATUS("Checking %s %c (ua %c) %5.2e %5.2e %5.2e %5.2f %5.2f %5.2f\n", + str_lattice(latt), cen, ua, a, b, c, al, be, ga); + + cell = cell_new_from_parameters(a, b, c, + deg2rad(al), deg2rad(be), deg2rad(ga)); + cell_set_lattice_type(cell, latt); + cell_set_centering(cell, cen); + cell_set_unique_axis(cell, ua); + + if ( check_cell(cell, "Input") ) fail = 1; + //cell_print(cell); + n = uncenter_cell(cell, &t); + if ( n != NULL ) { + if ( check_cell(n, "Output") ) fail = 1; + } else { + fail = 1; + } + + STATUS("Transformation was:\n"); + cell_transformation_print(t); - cell_set_parameters(cell, 10e-10, 20e-10, 30e-10, - deg2rad(90.0), deg2rad(90.0), deg2rad(100.0)); - cell_set_centering(cell, 'C'); - cell_set_lattice_type(cell, L_MONOCLINIC); - cell_set_unique_axis(cell, 'c'); - if ( check_cell(cell) ) fail = 1; - n = uncenter_cell(cell); - if ( check_cell(n) ) fail = 1; + if ( fail ) ERROR("\n"); return fail; } @@ -86,28 +103,100 @@ int main(int argc, char *argv[]) int fail = 0; /* Triclinic P */ - check_centering(10e-10, 20e-10, 30e-10, 100.0, 120.0, 140.0, - L_TRICLINIC, 'P', '*', - ); + fail += check_centering(50e-10, 55e-10, 70e-10, 67.0, 70.0, 77.0, + L_TRICLINIC, 'P', '*'); + /* Monoclinic P */ + fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0, + L_MONOCLINIC, 'P', 'a'); + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0, + L_MONOCLINIC, 'P', 'b'); + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0, + L_MONOCLINIC, 'P', 'c'); + /* Monoclinic A */ + fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0, + L_MONOCLINIC, 'A', 'a'); + /* Monoclinic B */ + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0, + L_MONOCLINIC, 'B', 'b'); + /* Monoclinic C */ + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0, + L_MONOCLINIC, 'C', 'c'); + /* Orthorhombic P */ + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_ORTHORHOMBIC, 'P', '*'); + /* Orthorhombic A */ + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_ORTHORHOMBIC, 'A', '*'); + /* Orthorhombic B */ + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_ORTHORHOMBIC, 'B', '*'); + /* Orthorhombic C */ + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_ORTHORHOMBIC, 'C', '*'); + /* Orthorhombic I */ + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_ORTHORHOMBIC, 'I', '*'); + /* Orthorhombic F */ + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_ORTHORHOMBIC, 'F', '*'); + /* Tetragonal P */ + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_TETRAGONAL, 'P', 'a'); + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_TETRAGONAL, 'P', 'b'); + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_TETRAGONAL, 'P', 'c'); + /* Tetragonal I */ + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_TETRAGONAL, 'I', 'a'); + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_TETRAGONAL, 'I', 'b'); + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, + L_TETRAGONAL, 'I', 'c'); + /* Rhombohedral R */ + fail += check_centering(10e-10, 10e-10, 10e-10, 60.0, 60.0, 60.0, + L_RHOMBOHEDRAL, 'R', '*'); + /* Hexagonal P */ + fail += check_centering(10e-10, 20e-10, 30e-10, 120.0, 90.0, 90.0, + L_HEXAGONAL, 'P', 'a'); + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 120.0, 90.0, + L_HEXAGONAL, 'P', 'b'); + fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 120.0, + L_HEXAGONAL, 'P', 'c'); + /* Hexagonal H (PDB-speak for rhombohedral) */ + fail += check_centering(40e-10, 20e-10, 20e-10, 120.0, 90.0, 90.0, + L_HEXAGONAL, 'H', 'a'); + fail += check_centering(20e-10, 40e-10, 20e-10, 90.0, 120.0, 90.0, + L_HEXAGONAL, 'H', 'b'); + fail += check_centering(20e-10, 20e-10, 40e-10, 90.0, 90.0, 120.0, + L_HEXAGONAL, 'H', 'c'); + /* Cubic P */ + fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0, + L_CUBIC, 'P', '*'); + /* Cubic I */ - /* Cubic F */ + fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0, + L_CUBIC, 'I', '*'); + /* Cubic F */ + fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0, + L_CUBIC, 'F', '*'); return fail; } -- cgit v1.2.3 From efd2a45e256d3d16d2f93ed24dedb8ecb9ca53ef Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 4 Sep 2012 18:33:24 +0200 Subject: WIP --- libcrystfel/src/cell-utils.c | 100 +++++++++++++++++-------------------------- libcrystfel/src/symmetry.c | 1 + 2 files changed, 41 insertions(+), 60 deletions(-) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 776526aa..b5d8f017 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -263,65 +263,61 @@ int bravais_lattice(UnitCell *cell) } +static UnitCellTransformation *tfn_identity() +{ +} + + +static void *tfn_combine(UnitCellTransformation *t, + double *na, double *nb, double *nc) +{ +} + + +static double *v(double a, double b, double c) +{ + double *vec = malloc(3*sizeof(double)); + if ( vec == NULL ) return NULL; + vec[0] = a; vec[1] = b; vec[2] = c; + return vec; +} + + static UnitCellTransformation *uncentering_transformation(UnitCell *in) { UnitCellTransformation *t; - 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; - 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); + t = tfn_identity(); + if ( t == NULL ) return NULL; + if ( (ua == 'a') || (cen == 'A') ) { - double tx, ty, tz; - tx = ax; ty = ay; tz = az; - ax = cx; ay = cy; az = cz; - cx = -tx; cy = -ty; cz = -tz; + tfn_combine(t, v(0,0,1), v(0,1,0), v(-1,0,0)); if ( cen == 'A' ) cen = 'C'; - ua = 'c'; } if ( (ua == 'b') || (cen == 'B') ) { - double tx, ty, tz; - tx = bx; ty = by; tz = bz; - bx = cx; by = cy; bz = cz; - cx = -tx; cy = -ty; cz = -tz; + tfn_combine(t, v(1,0,0), v(0,0,1), v(0,-1,0)); if ( cen == 'B' ) cen = 'C'; - ua = 'c'; } 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; + tfn_combine(t, v(-H,H,H), v(H,-H,H), v(H,H,-H)); + /* FIXME: How to handle the change of lattice type? */ if ( lt == L_CUBIC ) { cell_set_lattice_type(cell, L_RHOMBOHEDRAL); cell_set_centering(cell, 'R'); @@ -333,15 +329,7 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in) 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; + tfn_combine(t, v(0,H,H), v(H,0,H), v(H,H,0)); if ( lt == L_CUBIC ) { cell_set_lattice_type(cell, L_RHOMBOHEDRAL); cell_set_centering(cell, 'R'); @@ -354,30 +342,13 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in) 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; + tfn_combine(t, v(H,H,0), v(-H,H,0), v(0,0,1)); 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; + tfn_combine(t, v(TT,OT,OT), v(-OT,OT,OT), v(-OT,-TT,OT)); assert(lt == L_HEXAGONAL); cell_set_lattice_type(cell, L_RHOMBOHEDRAL); cell_set_centering(cell, 'R'); @@ -389,7 +360,17 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in) } + if ( ua == 'a' ) { + tfn_combine(t, v(0,0,-1), v(0,1,0), v(1,0,0)); + if ( cen == 'C' ) cen = 'A'; + } + if ( ua == 'b' ) { + tfn_combine(t, v(1,0,0), v(0,0,-1), v(0,1,0)); + if ( cen == 'C' ) cen = 'B'; + } + + return t; } @@ -408,7 +389,6 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in) */ UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t) { - UnitCell *cell; UnitCellTransformation *tt; if ( !bravais_lattice(in) ) { diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 68b39bcb..52a9aae6 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -256,6 +256,7 @@ int num_equivs(const SymOpList *ops, const SymOpMask *m) static signed int *v(signed int h, signed int k, signed int i, signed int l) { signed int *vec = malloc(3*sizeof(signed int)); + if ( vec == NULL ) return NULL; /* Convert back to 3-index form now */ vec[0] = h-i; vec[1] = k-i; vec[2] = l; return vec; -- cgit v1.2.3 From 42e1753bcfd692d5aaf6d7551c6e4635f65c683f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 6 Sep 2012 10:38:12 +0200 Subject: WIP --- libcrystfel/src/cell-utils.c | 20 -------------------- libcrystfel/src/cell.c | 40 +++++++++++++++++++++++++++++++++++++++- 2 files changed, 39 insertions(+), 21 deletions(-) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index b5d8f017..f0c6328f 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -263,26 +263,6 @@ int bravais_lattice(UnitCell *cell) } -static UnitCellTransformation *tfn_identity() -{ -} - - -static void *tfn_combine(UnitCellTransformation *t, - double *na, double *nb, double *nc) -{ -} - - -static double *v(double a, double b, double c) -{ - double *vec = malloc(3*sizeof(double)); - if ( vec == NULL ) return NULL; - vec[0] = a; vec[1] = b; vec[2] = c; - return vec; -} - - static UnitCellTransformation *uncentering_transformation(UnitCell *in) { UnitCellTransformation *t; diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index e2617dfc..53b0436d 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -634,7 +634,7 @@ UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t) /** - * cell_transform: + * cell_transform_inverse: * @cell: A %UnitCell. * @t: A %UnitCellTransformation. * @@ -649,6 +649,44 @@ UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t) } +/** + * tfn_identity: + * + * Returns: A %UnitCellTransformation corresponding to an identity operation. + * + */ +static UnitCellTransformation *tfn_identity() +{ +} + + +/** + * tfn_combine: + * @t: A %UnitCellTransformation + * @na: Pointer to three doubles representing naa, nab, nac + * @nb: Pointer to three doubles representing nba, nbb, nbc + * @nc: Pointer to three doubles representing nca, ncb, ncc + * + * Updates @t such that it represents its previous transformation followed by + * a new transformation, corresponding to letting a = naa*a + nab*b + nac*c. + * Likewise, a = nba*a + nbb*b + nbc*c and c = nca*a + ncb*b + ncc*c. + * + */ +static void tfn_combine(UnitCellTransformation *t, + double *na, double *nb, double *nc) +{ +} + + +static double *v(double a, double b, double c) +{ + double *vec = malloc(3*sizeof(double)); + if ( vec == NULL ) return NULL; + vec[0] = a; vec[1] = b; vec[2] = c; + return vec; +} + + void cell_transformation_print(UnitCellTransformation *t) { -- cgit v1.2.3 From 1363184f920a18fc560ad028a29ffcaa3148d93e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 7 Sep 2012 15:45:08 +0200 Subject: WIP --- libcrystfel/src/cell-utils.c | 101 ++++++++++++++++++++++++++++++++----------- libcrystfel/src/cell.c | 15 ++++--- libcrystfel/src/cell.h | 6 ++- tests/centering_check.c | 2 +- 4 files changed, 92 insertions(+), 32 deletions(-) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index f0c6328f..6c711642 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -263,7 +263,9 @@ int bravais_lattice(UnitCell *cell) } -static UnitCellTransformation *uncentering_transformation(UnitCell *in) +static UnitCellTransformation *uncentering_transformation(UnitCell *in, + char *new_centering, + LatticeType *new_latt) { UnitCellTransformation *t; const double OT = 1.0/3.0; @@ -280,12 +282,16 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in) if ( t == NULL ) return NULL; if ( (ua == 'a') || (cen == 'A') ) { - tfn_combine(t, v(0,0,1), v(0,1,0), v(-1,0,0)); + tfn_combine(t, tfn_vector(0,0,1), + tfn_vector(0,1,0), + tfn_vector(-1,0,0)); if ( cen == 'A' ) cen = 'C'; } if ( (ua == 'b') || (cen == 'B') ) { - tfn_combine(t, v(1,0,0), v(0,0,1), v(0,-1,0)); + tfn_combine(t, tfn_vector(1,0,0), + tfn_vector(0,0,1), + tfn_vector(0,-1,0)); if ( cen == 'B' ) cen = 'C'; } @@ -296,42 +302,50 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in) break; case 'I' : - tfn_combine(t, v(-H,H,H), v(H,-H,H), v(H,H,-H)); + tfn_combine(t, tfn_vector(-H,H,H), + tfn_vector(H,-H,H), + tfn_vector(H,H,-H)); /* FIXME: How to handle the change of lattice type? */ if ( lt == L_CUBIC ) { - cell_set_lattice_type(cell, L_RHOMBOHEDRAL); - cell_set_centering(cell, 'R'); + *new_latt = L_RHOMBOHEDRAL; + *new_centering = 'R'; } else { /* Tetragonal or orthorhombic */ - cell_set_lattice_type(cell, L_TRICLINIC); - cell_set_centering(cell, 'P'); + *new_latt = L_TRICLINIC; + *new_centering = 'P'; } break; case 'F' : - tfn_combine(t, v(0,H,H), v(H,0,H), v(H,H,0)); + tfn_combine(t, tfn_vector(0,H,H), + tfn_vector(H,0,H), + tfn_vector(H,H,0)); if ( lt == L_CUBIC ) { - cell_set_lattice_type(cell, L_RHOMBOHEDRAL); - cell_set_centering(cell, 'R'); + *new_latt = L_RHOMBOHEDRAL; + *new_centering = 'R'; } else { assert(lt == L_ORTHORHOMBIC); - cell_set_lattice_type(cell, L_TRICLINIC); - cell_set_centering(cell, 'P'); + *new_latt = L_TRICLINIC; + *new_centering = 'P'; } break; case 'C' : - tfn_combine(t, v(H,H,0), v(-H,H,0), v(0,0,1)); - cell_set_lattice_type(cell, L_MONOCLINIC); - cell_set_centering(cell, 'P'); + tfn_combine(t, tfn_vector(H,H,0), + tfn_vector(-H,H,0), + tfn_vector(0,0,1)); + *new_latt = L_MONOCLINIC; + *new_centering = 'P'; break; case 'H' : - tfn_combine(t, v(TT,OT,OT), v(-OT,OT,OT), v(-OT,-TT,OT)); + tfn_combine(t, tfn_vector(TT,OT,OT), + tfn_vector(-OT,OT,OT), + tfn_vector(-OT,-TT,OT)); assert(lt == L_HEXAGONAL); - cell_set_lattice_type(cell, L_RHOMBOHEDRAL); - cell_set_centering(cell, 'R'); + *new_latt = L_RHOMBOHEDRAL; + *new_centering = 'R'; break; default : @@ -341,12 +355,16 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in) } if ( ua == 'a' ) { - tfn_combine(t, v(0,0,-1), v(0,1,0), v(1,0,0)); + tfn_combine(t, tfn_vector(0,0,-1), + tfn_vector(0,1,0), + tfn_vector(1,0,0)); if ( cen == 'C' ) cen = 'A'; } if ( ua == 'b' ) { - tfn_combine(t, v(1,0,0), v(0,0,-1), v(0,1,0)); + tfn_combine(t, tfn_vector(1,0,0), + tfn_vector(0,0,-1), + tfn_vector(0,1,0)); if ( cen == 'C' ) cen = 'B'; } @@ -370,18 +388,27 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in) UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t) { UnitCellTransformation *tt; + char new_centering; + LatticeType new_latt; + UnitCell *out; if ( !bravais_lattice(in) ) { ERROR("Cannot uncenter: not a Bravais lattice.\n"); return NULL; } - tt = uncentering_transformation(in); + tt = uncentering_transformation(in, &new_centering, &new_latt); if ( tt == NULL ) return NULL; if ( t != NULL ) *t = tt; - return cell_transform(in, tt); + out = cell_transform(in, tt); + if ( out == NULL ) return NULL; + + cell_set_lattice_type(out, new_latt); + cell_set_centering(out, new_centering); + + return out; } @@ -634,12 +661,15 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, /* Reverse the de-centering transformation */ new_cell_trans = cell_transform_inverse(new_cell, to_given_cell); cell_free(new_cell); + cell_set_lattice_type(new_cell, cell_get_lattice_type(template_in)); + cell_set_centering(new_cell, cell_get_centering(template_in)); + cell_set_unique_axis(new_cell, cell_get_unique_axis(template_in)); return new_cell_trans; } -UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template) +UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in) { double ax, ay, az; double bx, by, bz; @@ -654,6 +684,18 @@ UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template) int have_real_a; int have_real_b; int have_real_c; + UnitCell *cell; + UnitCell *template; + UnitCellTransformation *to_given_cell; + UnitCell *new_cell; + UnitCell *new_cell_trans; + + /* "Un-center" the template unit cell to make the comparison easier */ + template = uncenter_cell(template_in, &to_given_cell); + + /* The candidate cell is also uncentered, because it might be centered + * if it came from (e.g.) MOSFLM */ + cell = uncenter_cell(cell_in, NULL); /* Get the lengths to match */ if ( cell_get_cartesian(template, &ax, &ay, &az, @@ -731,7 +773,16 @@ UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template) real_c.w = -real_c.w; } - return cell_new_from_direct_axes(real_a, real_b, real_c); + new_cell = cell_new_from_direct_axes(real_a, real_b, real_c); + + /* Reverse the de-centering transformation */ + new_cell_trans = cell_transform_inverse(new_cell, to_given_cell); + cell_free(new_cell); + cell_set_lattice_type(new_cell, cell_get_lattice_type(template_in)); + cell_set_centering(new_cell, cell_get_centering(template_in)); + cell_set_unique_axis(new_cell, cell_get_unique_axis(template_in)); + + return new_cell_trans; } diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 53b0436d..ca0fb9cd 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -601,6 +601,12 @@ const char *cell_rep(UnitCell *cell) } +struct _unitcelltransformation +{ + +}; + + static UnitCellTransformation *inverse_transformation(UnitCellTransformation *t) { /* FIXME: Implementation */ @@ -655,7 +661,7 @@ UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t) * Returns: A %UnitCellTransformation corresponding to an identity operation. * */ -static UnitCellTransformation *tfn_identity() +UnitCellTransformation *tfn_identity() { } @@ -672,13 +678,12 @@ static UnitCellTransformation *tfn_identity() * Likewise, a = nba*a + nbb*b + nbc*c and c = nca*a + ncb*b + ncc*c. * */ -static void tfn_combine(UnitCellTransformation *t, - double *na, double *nb, double *nc) +void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc) { } -static double *v(double a, double b, double c) +double *tfn_vector(double a, double b, double c) { double *vec = malloc(3*sizeof(double)); if ( vec == NULL ) return NULL; @@ -687,7 +692,7 @@ static double *v(double a, double b, double c) } -void cell_transformation_print(UnitCellTransformation *t) +void tfn_print(UnitCellTransformation *t) { } diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index 625884e1..e5c98ace 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -140,6 +140,10 @@ extern UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t); extern UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t); -extern void cell_transformation_print(UnitCellTransformation *t); +extern UnitCellTransformation *tfn_identity(void); +extern void tfn_combine(UnitCellTransformation *t, + double *na, double *nb, double *nc); +extern void tfn_print(UnitCellTransformation *t); +extern double *tfn_vector(double a, double b, double c); #endif /* CELL_H */ diff --git a/tests/centering_check.c b/tests/centering_check.c index 96326d21..bfee2806 100644 --- a/tests/centering_check.c +++ b/tests/centering_check.c @@ -89,7 +89,7 @@ static int check_centering(double a, double b, double c, } STATUS("Transformation was:\n"); - cell_transformation_print(t); + tfn_print(t); if ( fail ) ERROR("\n"); -- cgit v1.2.3 From ead022103e361778e150cce750db68b452f5753c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 11 Sep 2012 19:35:28 +0200 Subject: Implementation of UnitCellTransformation --- libcrystfel/src/cell.c | 158 +++++++++++++++++++++++++++++++++++++++++++++++-- libcrystfel/src/cell.h | 1 + 2 files changed, 153 insertions(+), 6 deletions(-) diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index ca0fb9cd..beaa341a 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -603,14 +603,46 @@ const char *cell_rep(UnitCell *cell) struct _unitcelltransformation { - + gsl_matrix *m; }; -static UnitCellTransformation *inverse_transformation(UnitCellTransformation *t) +UnitCellTransformation *tfn_inverse(UnitCellTransformation *t) { - /* FIXME: Implementation */ - return NULL; + int s; + gsl_matrix *inv; + gsl_permutation *perm; + UnitCellTransformation *out; + + out = tfn_identity(); + if ( out == NULL ) return NULL; + + perm = gsl_permutation_alloc(t->m->size1); + if ( perm == NULL ) { + ERROR("Couldn't allocate permutation\n"); + return NULL; + } + inv = gsl_matrix_alloc(t->m->size1, t->m->size2); + if ( inv == NULL ) { + ERROR("Couldn't allocate inverse\n"); + gsl_permutation_free(perm); + return NULL; + } + if ( gsl_linalg_LU_decomp(t->m, perm, &s) ) { + ERROR("Couldn't decompose matrix\n"); + gsl_permutation_free(perm); + return NULL; + } + if ( gsl_linalg_LU_invert(t->m, perm, inv) ) { + ERROR("Couldn't invert matrix\n"); + gsl_permutation_free(perm); + return NULL; + } + gsl_permutation_free(perm); + + gsl_matrix_free(out->m); + out->m = inv; + return out; } @@ -627,13 +659,52 @@ static UnitCellTransformation *inverse_transformation(UnitCellTransformation *t) UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t) { UnitCell *out; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + gsl_matrix *m; + gsl_matrix *a; if ( t == NULL ) return NULL; out = cell_new(); if ( out == NULL ) return NULL; - /* FIXME: Implementation */ + cell_get_cartesian(cell, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + + m = gsl_matrix_alloc(3,3); + a = gsl_matrix_alloc(3,3); + if ( (m == NULL) || (a == NULL) ) { + cell_free(out); + return NULL; + } + + gsl_matrix_set(m, 0, 0, ax); + gsl_matrix_set(m, 0, 1, ay); + gsl_matrix_set(m, 0, 2, az); + gsl_matrix_set(m, 1, 0, bx); + gsl_matrix_set(m, 1, 1, by); + gsl_matrix_set(m, 1, 2, bz); + gsl_matrix_set(m, 2, 0, cx); + gsl_matrix_set(m, 2, 1, cy); + gsl_matrix_set(m, 2, 2, cz); + + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, t->m, m, 1.0, a); + + cell_set_cartesian(out, gsl_matrix_get(a, 0, 0), + gsl_matrix_get(a, 0, 1), + gsl_matrix_get(a, 0, 2), + gsl_matrix_get(a, 1, 0), + gsl_matrix_get(a, 1, 1), + gsl_matrix_get(a, 1, 2), + gsl_matrix_get(a, 2, 0), + gsl_matrix_get(a, 2, 1), + gsl_matrix_get(a, 2, 2)); + + gsl_matrix_free(a); + gsl_matrix_free(m); return out; } @@ -651,7 +722,7 @@ UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t) */ UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t) { - return cell_transform(cell, inverse_transformation(t)); + return cell_transform(cell, tfn_inverse(t)); } @@ -663,6 +734,20 @@ UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t) */ UnitCellTransformation *tfn_identity() { + UnitCellTransformation *tfn; + + tfn = calloc(1, sizeof(UnitCellTransformation)); + if ( tfn == NULL ) return NULL; + + tfn->m = gsl_matrix_alloc(3, 3); + if ( tfn->m == NULL ) { + free(tfn); + return NULL; + } + + gsl_matrix_set_identity(tfn->m); + + return tfn; } @@ -680,6 +765,33 @@ UnitCellTransformation *tfn_identity() */ void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc) { + gsl_matrix *a; + gsl_matrix *n; + + n = gsl_matrix_alloc(3, 3); + a = gsl_matrix_alloc(3, 3); + if ( (n == NULL) || (a == NULL) ) { + return; + } + gsl_matrix_set(n, 0, 0, na[0]); + gsl_matrix_set(n, 0, 1, na[1]); + gsl_matrix_set(n, 0, 2, na[2]); + gsl_matrix_set(n, 1, 0, nb[0]); + gsl_matrix_set(n, 1, 1, nb[1]); + gsl_matrix_set(n, 1, 2, nb[2]); + gsl_matrix_set(n, 2, 0, nc[0]); + gsl_matrix_set(n, 2, 1, nc[1]); + gsl_matrix_set(n, 2, 2, nc[2]); + + free(na); + free(nb); + free(nc); + + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, n, t->m, 1.0, a); + + gsl_matrix_free(t->m); + t->m = a; + gsl_matrix_free(n); } @@ -694,5 +806,39 @@ double *tfn_vector(double a, double b, double c) void tfn_print(UnitCellTransformation *t) { + gsl_permutation *perm; + gsl_matrix *lu; + int s; + + STATUS("New a = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 0, 0), + gsl_matrix_get(t->m, 0, 1), + gsl_matrix_get(t->m, 0, 2)); + STATUS("New b = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 1, 0), + gsl_matrix_get(t->m, 1, 1), + gsl_matrix_get(t->m, 1, 2)); + STATUS("New c = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 2, 0), + gsl_matrix_get(t->m, 2, 1), + gsl_matrix_get(t->m, 2, 2)); + lu = gsl_matrix_alloc(3, 3); + if ( lu == NULL ) { + ERROR("Couldn't allocate LU decomposition.\n"); + return; + } + + gsl_matrix_memcpy(lu, t->m); + + perm = gsl_permutation_alloc(t->m->size1); + if ( perm == NULL ) { + ERROR("Couldn't allocate permutation.\n"); + gsl_matrix_free(lu); + return; + } + if ( gsl_linalg_LU_decomp(lu, perm, &s) ) { + ERROR("LU decomposition failed.\n"); + gsl_permutation_free(perm); + gsl_matrix_free(lu); + return; + } + STATUS("Transformation determinant = %.2f\n", gsl_linalg_LU_det(lu, s)); } diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index e5c98ace..ac1bd555 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -144,6 +144,7 @@ extern UnitCellTransformation *tfn_identity(void); extern void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc); extern void tfn_print(UnitCellTransformation *t); +extern UnitCellTransformation *tfn_inverse(UnitCellTransformation *t); extern double *tfn_vector(double a, double b, double c); #endif /* CELL_H */ -- cgit v1.2.3 From bf986af745b88ddd0c1309c3af03d06165a6c95f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 11 Sep 2012 19:35:38 +0200 Subject: Remove unnecessary includes --- libcrystfel/src/reax.c | 1 - src/hrs-scaling.c | 1 - src/post-refinement.c | 1 - 3 files changed, 3 deletions(-) diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c index 29ea7294..7529d12f 100644 --- a/libcrystfel/src/reax.c +++ b/libcrystfel/src/reax.c @@ -41,7 +41,6 @@ #include #include #include -#include #include "image.h" #include "utils.h" diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 3ae0d563..399d5f90 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -37,7 +37,6 @@ #include #include #include -#include #include "image.h" #include "peaks.h" diff --git a/src/post-refinement.c b/src/post-refinement.c index ac4de1bc..7410d931 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -37,7 +37,6 @@ #include #include #include -#include #include "image.h" #include "post-refinement.h" -- cgit v1.2.3 From 5911a4b627e1c6676c5522a8ecd61a2a834f544e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 12 Sep 2012 17:57:38 +0200 Subject: More work on transformations --- Makefile.am | 8 +++- doc/reference/CrystFEL-docs.sgml | 1 + doc/reference/CrystFEL-sections.txt | 16 ++++++++ libcrystfel/src/cell-utils.c | 13 ++++++ libcrystfel/src/cell.c | 81 ++++++++++++++++++++++++++++++++----- libcrystfel/src/cell.h | 1 + tests/centering_check.c | 3 +- tests/transformation_check.c | 75 ++++++++++++++++++++++++++++++++++ 8 files changed, 186 insertions(+), 12 deletions(-) create mode 100644 tests/transformation_check.c diff --git a/Makefile.am b/Makefile.am index 953a04df..f0adc640 100644 --- a/Makefile.am +++ b/Makefile.am @@ -8,12 +8,12 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \ noinst_PROGRAMS = tests/list_check tests/integration_check \ tests/pr_gradient_check tests/symmetry_check \ - tests/centering_check + tests/centering_check tests/transformation_check TESTS = tests/list_check tests/first_merge_check tests/second_merge_check \ tests/third_merge_check tests/fourth_merge_check \ tests/integration_check tests/pr_gradient_check \ - tests/symmetry_check tests/centering_check + tests/symmetry_check tests/centering_check tests/transformation_check EXTRA_DIST += tests/first_merge_check tests/second_merge_check \ tests/third_merge_check tests/fourth_merge_check @@ -87,6 +87,10 @@ tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c \ tests_centering_check_SOURCES = tests/centering_check.c libcrystfel/src/cell.c \ libcrystfel/src/cell-utils.c +tests_transformation_check_SOURCES = tests/transformation_check.c \ + libcrystfel/src/cell.c \ + libcrystfel/src/cell-utils.c + INCLUDES = -I$(top_srcdir)/libcrystfel/src -I$(top_srcdir)/data EXTRA_DIST += src/dw-hdfsee.h src/hdfsee.h src/render_hkl.h \ diff --git a/doc/reference/CrystFEL-docs.sgml b/doc/reference/CrystFEL-docs.sgml index 89b4832e..e1da1784 100644 --- a/doc/reference/CrystFEL-docs.sgml +++ b/doc/reference/CrystFEL-docs.sgml @@ -34,6 +34,7 @@ Unit cells + diff --git a/doc/reference/CrystFEL-sections.txt b/doc/reference/CrystFEL-sections.txt index 284740a1..3ccc1511 100644 --- a/doc/reference/CrystFEL-sections.txt +++ b/doc/reference/CrystFEL-sections.txt @@ -88,12 +88,28 @@ cell_set_pointgroup cell_set_reciprocal cell_set_spacegroup +tfn_combine +tfn_identity +tfn_inverse +tfn_print +tfn_vector +tfn_free + + +
+cell-utils cell_rotate rotate_cell cell_print resolution match_cell match_cell_ab +cell_is_sensible +validate_cell +uncenter_cell +bravais_lattice +right_handed +str_lattice
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 6c711642..9b65cae3 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -47,6 +47,19 @@ #include "image.h" +/** + * SECTION:cell-utils + * @short_description: Unit cell utilities + * @title: Unit cell utilities + * @section_id: + * @see_also: + * @include: "cell-utils.h" + * @Image: + * + * There are some utility functions associated with the core %UnitCell. + **/ + + /* Weighting factor of lengths relative to angles */ #define LWEIGHT (10.0e-9) diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index beaa341a..7a62f51c 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -607,33 +607,55 @@ struct _unitcelltransformation }; +/** + * tfn_inverse: + * @t: A %UnitCellTransformation. + * + * Calculates the inverse of @t. That is, if you apply cell_transform() to a + * %UnitCell using @t, and then apply cell_transform() to the result using + * tfn_inverse(@t) instead of @t, you will recover the same lattice vectors + * (but note that the lattice type, centering and unique axis information will + * be lost). + * + * Returns: The inverse of @t. + * + */ UnitCellTransformation *tfn_inverse(UnitCellTransformation *t) { int s; + gsl_matrix *m; gsl_matrix *inv; gsl_permutation *perm; UnitCellTransformation *out; + m = gsl_matrix_alloc(3, 3); + if ( m == NULL ) return NULL; + out = tfn_identity(); - if ( out == NULL ) return NULL; + if ( out == NULL ) { + gsl_matrix_free(m); + return NULL; + } - perm = gsl_permutation_alloc(t->m->size1); + gsl_matrix_memcpy(m, t->m); + + perm = gsl_permutation_alloc(m->size1); if ( perm == NULL ) { ERROR("Couldn't allocate permutation\n"); return NULL; } - inv = gsl_matrix_alloc(t->m->size1, t->m->size2); + inv = gsl_matrix_alloc(m->size1, m->size2); if ( inv == NULL ) { ERROR("Couldn't allocate inverse\n"); gsl_permutation_free(perm); return NULL; } - if ( gsl_linalg_LU_decomp(t->m, perm, &s) ) { + if ( gsl_linalg_LU_decomp(m, perm, &s) ) { ERROR("Couldn't decompose matrix\n"); gsl_permutation_free(perm); return NULL; } - if ( gsl_linalg_LU_invert(t->m, perm, inv) ) { + if ( gsl_linalg_LU_invert(m, perm, inv) ) { ERROR("Couldn't invert matrix\n"); gsl_permutation_free(perm); return NULL; @@ -641,6 +663,7 @@ UnitCellTransformation *tfn_inverse(UnitCellTransformation *t) gsl_permutation_free(perm); gsl_matrix_free(out->m); + gsl_matrix_free(m); out->m = inv; return out; } @@ -651,7 +674,8 @@ UnitCellTransformation *tfn_inverse(UnitCellTransformation *t) * @cell: A %UnitCell. * @t: A %UnitCellTransformation. * - * Applies @t to @cell. + * Applies @t to @cell. Note that the lattice type, centering and unique axis + * information will not be preserved. * * Returns: Transformed copy of @cell. * @@ -675,7 +699,7 @@ UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t) &cx, &cy, &cz); m = gsl_matrix_alloc(3,3); - a = gsl_matrix_alloc(3,3); + a = gsl_matrix_calloc(3,3); if ( (m == NULL) || (a == NULL) ) { cell_free(out); return NULL; @@ -722,7 +746,13 @@ UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t) */ UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t) { - return cell_transform(cell, tfn_inverse(t)); + UnitCellTransformation *inv; + UnitCell *out; + + inv = tfn_inverse(t); + out = cell_transform(cell, inv); + tfn_free(inv); + return out; } @@ -769,7 +799,7 @@ void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc) gsl_matrix *n; n = gsl_matrix_alloc(3, 3); - a = gsl_matrix_alloc(3, 3); + a = gsl_matrix_calloc(3, 3); if ( (n == NULL) || (a == NULL) ) { return; } @@ -795,6 +825,18 @@ void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc) } +/** + * tfn_vector: + * @a: Amount of "a" to include in new vector + * @b: Amount of "b" to include in new vector + * @c: Amount of "c" to include in new vector + * + * This is a convenience function to use when sending vectors to tfn_combine(): + * tfn_combine(tfn, tfn_vector(1,0,0), + * tfn_vector(0,2,0), + * tfn_vector(0,0,1)); + * + */ double *tfn_vector(double a, double b, double c) { double *vec = malloc(3*sizeof(double)); @@ -804,6 +846,13 @@ double *tfn_vector(double a, double b, double c) } +/** + * tfn_print: + * @t: A %UnitCellTransformation + * + * Prints information about @t to stderr. + * + */ void tfn_print(UnitCellTransformation *t) { gsl_permutation *perm; @@ -842,3 +891,17 @@ void tfn_print(UnitCellTransformation *t) STATUS("Transformation determinant = %.2f\n", gsl_linalg_LU_det(lu, s)); } + + +/** + * tfn_free: + * @t: A %UnitCellTransformation + * + * Frees all resources associated with @t. + * + */ +void tfn_free(UnitCellTransformation *t) +{ + gsl_matrix_free(t->m); + free(t); +} diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index ac1bd555..c7b8f8d6 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -146,5 +146,6 @@ extern void tfn_combine(UnitCellTransformation *t, extern void tfn_print(UnitCellTransformation *t); extern UnitCellTransformation *tfn_inverse(UnitCellTransformation *t); extern double *tfn_vector(double a, double b, double c); +extern void tfn_free(UnitCellTransformation *t); #endif /* CELL_H */ diff --git a/tests/centering_check.c b/tests/centering_check.c index bfee2806..eca67c85 100644 --- a/tests/centering_check.c +++ b/tests/centering_check.c @@ -70,7 +70,8 @@ static int check_centering(double a, double b, double c, UnitCellTransformation *t; int fail = 0; - STATUS("Checking %s %c (ua %c) %5.2e %5.2e %5.2e %5.2f %5.2f %5.2f\n", + STATUS(" ---------------> " + "Checking %s %c (ua %c) %5.2e %5.2e %5.2e %5.2f %5.2f %5.2f\n", str_lattice(latt), cen, ua, a, b, c, al, be, ga); cell = cell_new_from_parameters(a, b, c, diff --git a/tests/transformation_check.c b/tests/transformation_check.c new file mode 100644 index 00000000..771e44fe --- /dev/null +++ b/tests/transformation_check.c @@ -0,0 +1,75 @@ +/* + * transformation_check.c + * + * Check that unit cell transformations work + * + * Copyright © 2012 Thomas White + * + * 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 + + +int main(int argc, char *argv[]) +{ + int fail = 0; + UnitCell *cell, *cnew, *cback; + UnitCellTransformation *tfn, *inv; + + cell = cell_new_from_parameters(50e-10, 55e-10, 70e-10, + deg2rad(67.0), + deg2rad(70.0), + deg2rad(77.0)); + if ( cell == NULL ) return 1; + + tfn = tfn_identity(); + if ( tfn == NULL ) return 1; + + tfn_combine(tfn, tfn_vector(0,1,0), + tfn_vector(1,0,0), + tfn_vector(0,0,1)); + + + cell_print(cell); + tfn_print(tfn); + + cnew = cell_transform(cell, tfn); + cell_print(cnew); + + cback = cell_transform_inverse(cnew, tfn); + inv = tfn_inverse(tfn); + tfn_print(inv); + cell_print(cback); + + cell_get_cartesian(cell, &ax1, &ay1, &az1, + &by1, &by1, &bz1, + &cx1, &cy1, &cz1); + cell_get_cartesian(cback, &ax, &ay, &az, &by, &by, &bz, &cx, &cy, &cz); + + return fail; +} -- cgit v1.2.3 From 5b20bb5f7c5815ffa1aad3080269c519e0b8283f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 12 Sep 2012 17:58:15 +0200 Subject: Ignore new test --- tests/.gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/.gitignore b/tests/.gitignore index 16acffc3..81fd558f 100644 --- a/tests/.gitignore +++ b/tests/.gitignore @@ -6,4 +6,5 @@ integration_check pr_gradient_check symmetry_check centering_check +transformation_check .dirstamp -- cgit v1.2.3 From a670d325edde9b4580ef31afaad988456660c4ff Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 26 Sep 2012 12:18:19 +0200 Subject: Finish transformation_check --- tests/transformation_check.c | 94 +++++++++++++++++++++++++++++++++++--------- 1 file changed, 75 insertions(+), 19 deletions(-) diff --git a/tests/transformation_check.c b/tests/transformation_check.c index 771e44fe..80324795 100644 --- a/tests/transformation_check.c +++ b/tests/transformation_check.c @@ -35,41 +35,97 @@ #include +static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn) +{ + UnitCell *cnew, *cback; + UnitCellTransformation *inv; + double a[9], b[9]; + int i; + int fail = 0; + + cnew = cell_transform(cell, tfn); + + cback = cell_transform_inverse(cnew, tfn); + + cell_get_cartesian(cell, &a[0], &a[1], &a[2], + &a[3], &a[4], &a[5], + &a[6], &a[7], &a[8]); + cell_get_cartesian(cback, &b[0], &b[1], &b[2], + &b[3], &b[4], &b[5], + &b[6], &b[7], &b[8]); + for ( i=0; i<9; i++ ) { + if ( !within_tolerance(a[i], b[i], 0.1) ) { + fail = 1; + STATUS("%e %e\n", a[i], b[i]); + } + } + + if ( fail ) { + ERROR("Original cell not recovered after transformation:\n"); + cell_print(cell); + tfn_print(tfn); + inv = tfn_inverse(tfn); + tfn_print(inv); + cell_print(cback); + } + + return fail; +} + + int main(int argc, char *argv[]) { int fail = 0; - UnitCell *cell, *cnew, *cback; - UnitCellTransformation *tfn, *inv; + UnitCell *cell, *cref; + UnitCellTransformation *tfn; - cell = cell_new_from_parameters(50e-10, 55e-10, 70e-10, + cref = cell_new_from_parameters(50e-10, 55e-10, 70e-10, deg2rad(67.0), deg2rad(70.0), deg2rad(77.0)); + if ( cref == NULL ) return 1; + + cell = cell_rotate(cref, random_quaternion()); if ( cell == NULL ) return 1; + cell_free(cref); + /* Permutation of axes */ tfn = tfn_identity(); if ( tfn == NULL ) return 1; - tfn_combine(tfn, tfn_vector(0,1,0), - tfn_vector(1,0,0), - tfn_vector(0,0,1)); + tfn_vector(0,0,1), + tfn_vector(1,0,0)); + fail += check_transformation(cell, tfn); + tfn_free(tfn); + /* Doubling of cell in one direction */ + tfn = tfn_identity(); + if ( tfn == NULL ) return 1; + tfn_combine(tfn, tfn_vector(2,0,0), + tfn_vector(0,1,0), + tfn_vector(0,0,1)); + fail += check_transformation(cell, tfn); + tfn_free(tfn); - cell_print(cell); - tfn_print(tfn); + /* Diagonal supercell */ + tfn = tfn_identity(); + if ( tfn == NULL ) return 1; + tfn_combine(tfn, tfn_vector(1,1,0), + tfn_vector(0,1,0), + tfn_vector(0,0,1)); + fail += check_transformation(cell, tfn); + tfn_free(tfn); - cnew = cell_transform(cell, tfn); - cell_print(cnew); + /* Crazy */ + tfn = tfn_identity(); + if ( tfn == NULL ) return 1; + tfn_combine(tfn, tfn_vector(1,1,0), + tfn_vector(0,1,0), + tfn_vector(0,1,1)); + fail += check_transformation(cell, tfn); + tfn_free(tfn); - cback = cell_transform_inverse(cnew, tfn); - inv = tfn_inverse(tfn); - tfn_print(inv); - cell_print(cback); - - cell_get_cartesian(cell, &ax1, &ay1, &az1, - &by1, &by1, &bz1, - &cx1, &cy1, &cz1); - cell_get_cartesian(cback, &ax, &ay, &az, &by, &by, &bz, &cx, &cy, &cz); + cell_free(cell); return fail; } -- cgit v1.2.3 From 7b31c8f41732b6ca5402afd8fe6cb4d6e185bc48 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 26 Sep 2012 15:32:32 +0200 Subject: Formatting --- libcrystfel/src/cell-utils.c | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 9b65cae3..2f966466 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -315,10 +315,9 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in, break; case 'I' : - tfn_combine(t, tfn_vector(-H,H,H), - tfn_vector(H,-H,H), - tfn_vector(H,H,-H)); - /* FIXME: How to handle the change of lattice type? */ + tfn_combine(t, tfn_vector(-H,H,H), + tfn_vector(H,-H,H), + tfn_vector(H,H,-H)); if ( lt == L_CUBIC ) { *new_latt = L_RHOMBOHEDRAL; *new_centering = 'R'; @@ -330,9 +329,9 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in, break; case 'F' : - tfn_combine(t, tfn_vector(0,H,H), - tfn_vector(H,0,H), - tfn_vector(H,H,0)); + tfn_combine(t, tfn_vector(0,H,H), + tfn_vector(H,0,H), + tfn_vector(H,H,0)); if ( lt == L_CUBIC ) { *new_latt = L_RHOMBOHEDRAL; *new_centering = 'R'; -- cgit v1.2.3 From 487f375b726fd555ff0c7eab40499201fc2e1d7d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 27 Sep 2012 12:28:57 +0200 Subject: Finishing tweaks for uncentering stuff Today, I have mostly been having my life made difficult by the PDB's invention of "H centering". --- libcrystfel/src/cell-utils.c | 128 ++++++++++++++++++++++++++-------- libcrystfel/src/cell-utils.h | 5 +- libcrystfel/src/cell.c | 11 +-- tests/centering_check.c | 159 ++++++++++++++++++++++++++++++++++--------- tests/transformation_check.c | 45 ++++++++++++ 5 files changed, 280 insertions(+), 68 deletions(-) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 2f966466..f854efe0 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -171,21 +171,19 @@ void cell_print(UnitCell *cell) double a, b, c, alpha, beta, gamma; double ax, ay, az, bx, by, bz, cx, cy, cz; LatticeType lt; + char cen; lt = cell_get_lattice_type(cell); + cen = cell_get_centering(cell); - STATUS("%s %c", str_lattice(lt), cell_get_centering(cell)); + STATUS("%s %c", str_lattice(lt), cen); - switch ( lt ) + if ( (lt==L_MONOCLINIC) || (lt==L_TETRAGONAL) || ( lt==L_HEXAGONAL) + || ( (lt==L_ORTHORHOMBIC) && (cen=='A') ) + || ( (lt==L_ORTHORHOMBIC) && (cen=='B') ) + || ( (lt==L_ORTHORHOMBIC) && (cen=='C') ) ) { - 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) ) { @@ -263,6 +261,11 @@ int bravais_lattice(UnitCell *cell) return 0; case 'H' : + /* "Hexagonal H" is not a Bravais lattice, but rather something + * invented by the PDB to make life difficult for programmers. + * Accepting it as Bravais seems to be the least painful way to + * handle it correctly. Yuk. */ + if ( ua != 'c' ) return 0; if ( lattice == L_HEXAGONAL ) return 1; return 0; @@ -294,24 +297,28 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in, t = tfn_identity(); if ( t == NULL ) return NULL; - if ( (ua == 'a') || (cen == 'A') ) { + if ( ua == 'a' ) { tfn_combine(t, tfn_vector(0,0,1), tfn_vector(0,1,0), tfn_vector(-1,0,0)); - if ( cen == 'A' ) cen = 'C'; } - if ( (ua == 'b') || (cen == 'B') ) { + if ( ua == 'b' ) { tfn_combine(t, tfn_vector(1,0,0), tfn_vector(0,0,1), tfn_vector(0,-1,0)); - if ( cen == 'B' ) cen = 'C'; } switch ( cen ) { case 'P' : + *new_latt = lt; + *new_centering = 'P'; + break; + case 'R' : + *new_latt = L_RHOMBOHEDRAL; + *new_centering = 'R'; break; case 'I' : @@ -335,7 +342,6 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in, if ( lt == L_CUBIC ) { *new_latt = L_RHOMBOHEDRAL; *new_centering = 'R'; - } else { assert(lt == L_ORTHORHOMBIC); *new_latt = L_TRICLINIC; @@ -343,6 +349,8 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in, } break; + case 'A' : + case 'B' : case 'C' : tfn_combine(t, tfn_vector(H,H,0), tfn_vector(-H,H,0), @@ -352,6 +360,7 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in, break; case 'H' : + /* Obverse setting */ tfn_combine(t, tfn_vector(TT,OT,OT), tfn_vector(-OT,OT,OT), tfn_vector(-OT,-TT,OT)); @@ -366,18 +375,20 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in, } - if ( ua == 'a' ) { - tfn_combine(t, tfn_vector(0,0,-1), - tfn_vector(0,1,0), - tfn_vector(1,0,0)); - if ( cen == 'C' ) cen = 'A'; - } + /* Reverse the axis permutation, but only if this was not an H->R + * transformation */ + if ( !((cen=='H') && (*new_latt == L_RHOMBOHEDRAL)) ) { + if ( ua == 'a' ) { + tfn_combine(t, tfn_vector(0,0,-1), + tfn_vector(0,1,0), + tfn_vector(1,0,0)); + } - if ( ua == 'b' ) { - tfn_combine(t, tfn_vector(1,0,0), - tfn_vector(0,0,-1), - tfn_vector(0,1,0)); - if ( cen == 'C' ) cen = 'B'; + if ( ua == 'b' ) { + tfn_combine(t, tfn_vector(1,0,0), + tfn_vector(0,0,-1), + tfn_vector(0,1,0)); + } } return t; @@ -406,6 +417,7 @@ UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t) if ( !bravais_lattice(in) ) { ERROR("Cannot uncenter: not a Bravais lattice.\n"); + cell_print(in); return NULL; } @@ -481,11 +493,11 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, float angtol = deg2rad(tols[3]); UnitCell *cell; UnitCell *template; - UnitCellTransformation *to_given_cell; + UnitCellTransformation *uncentering; UnitCell *new_cell_trans; /* "Un-center" the template unit cell to make the comparison easier */ - template = uncenter_cell(template_in, &to_given_cell); + template = uncenter_cell(template_in, &uncentering); /* The candidate cell is also uncentered, because it might be centered * if it came from (e.g.) MOSFLM */ @@ -671,7 +683,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, free(cand[2]); /* Reverse the de-centering transformation */ - new_cell_trans = cell_transform_inverse(new_cell, to_given_cell); + new_cell_trans = cell_transform_inverse(new_cell, uncentering); cell_free(new_cell); cell_set_lattice_type(new_cell, cell_get_lattice_type(template_in)); cell_set_centering(new_cell, cell_get_centering(template_in)); @@ -1115,10 +1127,13 @@ int cell_is_sensible(UnitCell *cell) * lattice is a conventional Bravais lattice. * Warnings are printied if any of the checks are failed. * + * Returns: true if cell is invalid. + * */ -void validate_cell(UnitCell *cell) +int validate_cell(UnitCell *cell) { int err = 0; + char cen, ua; if ( !cell_is_sensible(cell) ) { ERROR("Warning: Unit cell parameters are not sensible.\n"); @@ -1136,5 +1151,58 @@ void validate_cell(UnitCell *cell) err = 1; } - if ( err ) cell_print(cell); + cen = cell_get_centering(cell); + ua = cell_get_unique_axis(cell); + if ( (cen == 'A') && (ua != 'a') ) { + ERROR("Warning: centering doesn't match unique axis.\n"); + err = 1; + } + if ( (cen == 'B') && (ua != 'b') ) { + ERROR("Warning: centering doesn't match unique axis.\n"); + err = 1; + } + if ( (cen == 'C') && (ua != 'c') ) { + ERROR("Warning: centering doesn't match unique axis.\n"); + err = 1; + } + + return err; +} + + +/** + * forbidden_reflection: + * @cell: A %UnitCell + * @h: h index to check + * @k: k index to check + * @l: l index to check + * + * Returns: true if this reflection is forbidden. + * + */ +int forbidden_reflection(UnitCell *cell, + signed int h, signed int k, signed int l) +{ + char cen; + + cen = cell_get_centering(cell); + + /* Reflection conditions here must match the transformation matrices + * in uncentering_transformation(). tests/centering_check verifies + * this (amongst other things). */ + + if ( cen == 'P' ) return 0; + if ( cen == 'R' ) return 0; + + if ( cen == 'A' ) return (k+l) % 2; + if ( cen == 'B' ) return (h+l) % 2; + if ( cen == 'C' ) return (h+k) % 2; + + if ( cen == 'I' ) return (h+k+l) % 2; + if ( cen == 'F' ) return ((h+k) % 2) || ((h+l) % 2) || ((k+l) % 2); + + /* Obverse setting */ + if ( cen == 'H' ) return (-h+k+l) % 3; + + return 0; } diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index 8acf2f85..f92ab22d 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -53,7 +53,7 @@ extern UnitCell *load_cell_from_pdb(const char *filename); extern int cell_is_sensible(UnitCell *cell); -extern void validate_cell(UnitCell *cell); +extern int validate_cell(UnitCell *cell); extern UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **tr); @@ -63,4 +63,7 @@ extern int right_handed(UnitCell *cell); extern const char *str_lattice(LatticeType l); +extern int forbidden_reflection(UnitCell *cell, + signed int h, signed int k, signed int l); + #endif /* CELL_UTILS_H */ diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 7a62f51c..e7c9dace 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -260,6 +260,9 @@ UnitCell *cell_new_from_cell(UnitCell *orig) cell_get_cartesian(orig, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); cell_set_cartesian(new, ax, ay, az, bx, by, bz, cx, cy, cz); cell_set_pointgroup(new, orig->pointgroup); + cell_set_lattice_type(new, orig->lattice_type); + cell_set_centering(new, orig->centering); + cell_set_unique_axis(new, orig->unique_axis); return new; } @@ -691,12 +694,12 @@ UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t) if ( t == NULL ) return NULL; - out = cell_new(); + out = cell_new_from_cell(cell); if ( out == NULL ) return NULL; - cell_get_cartesian(cell, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz); + cell_get_cartesian(out, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); m = gsl_matrix_alloc(3,3); a = gsl_matrix_calloc(3,3); diff --git a/tests/centering_check.c b/tests/centering_check.c index eca67c85..e17915d3 100644 --- a/tests/centering_check.c +++ b/tests/centering_check.c @@ -30,6 +30,7 @@ #include #include #include +#include #include #include @@ -39,24 +40,13 @@ static int check_cell(UnitCell *cell, const char *text) { int err = 0; - if ( !cell_is_sensible(cell) ) { - ERROR(" %s unit cell parameters are not sensible.\n", text); - err = 1; - } - - if ( !bravais_lattice(cell) ) { - ERROR(" %s unit cell is not a conventional Bravais" - " lattice.\n", text); - err = 1; - } + err = validate_cell(cell); - if ( !right_handed(cell) ) { - ERROR(" %s unit cell is not right handed.\n", text); - err = 1; + if ( err ) { + ERROR("%s cell:\n", text); + cell_print(cell); } - if ( err ) cell_print(cell); - return err; } @@ -65,40 +55,147 @@ static int check_centering(double a, double b, double c, double al, double be, double ga, LatticeType latt, char cen, char ua) { - UnitCell *cell; + UnitCell *cell, *cref; UnitCell *n; UnitCellTransformation *t; int fail = 0; + int i; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; STATUS(" ---------------> " "Checking %s %c (ua %c) %5.2e %5.2e %5.2e %5.2f %5.2f %5.2f\n", str_lattice(latt), cen, ua, a, b, c, al, be, ga); - cell = cell_new_from_parameters(a, b, c, + cref = cell_new_from_parameters(a, b, c, deg2rad(al), deg2rad(be), deg2rad(ga)); - cell_set_lattice_type(cell, latt); - cell_set_centering(cell, cen); - cell_set_unique_axis(cell, ua); + cell_set_lattice_type(cref, latt); + cell_set_centering(cref, cen); + cell_set_unique_axis(cref, ua); + + cell = cell_rotate(cref, random_quaternion()); + if ( cell == NULL ) return 1; + cell_free(cref); - if ( check_cell(cell, "Input") ) fail = 1; - //cell_print(cell); + check_cell(cell, "Input"); n = uncenter_cell(cell, &t); if ( n != NULL ) { + STATUS("Transformation was:\n"); + tfn_print(t); if ( check_cell(n, "Output") ) fail = 1; + if ( !fail ) cell_print(n); } else { fail = 1; } - STATUS("Transformation was:\n"); - tfn_print(t); + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + cell_get_cartesian(n, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + + fesetround(1); /* Round towards nearest */ + for ( i=0; i<100; i++ ) { + + signed int h, k, l; + double x, y, z; + double nh, nk, nl; + double dh, dk, dl; + int f = 0; + + do { + + h = flat_noise(0, 30); + k = flat_noise(0, 30); + l = flat_noise(0, 30); + + } while ( forbidden_reflection(cell, h, k, l) ); - if ( fail ) ERROR("\n"); + x = h*asx + k*bsx + l*csx; + y = h*asy + k*bsy + l*csy; + z = h*asz + k*bsz + l*csz; + + nh = x*ax + y*ay + z*az; + nk = x*bx + y*by + z*bz; + nl = x*cx + y*cy + z*cz; + + dh = nh - lrint(nh); + dk = nk - lrint(nk); + dl = nl - lrint(nl); + if ( fabs(dh) > 0.1 ) f++; + if ( fabs(dk) > 0.1 ) f++; + if ( fabs(dl) > 0.1 ) f++; + + if ( f ) { + STATUS("Centered %3i %3i %3i -> " + "Primitive %7.2f %7.2f %7.2f\n", + h, k, l, nh, nk, nl); + fail = 1; + } + + } + + cell_get_reciprocal(n, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + cell_get_cartesian(cell, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + + for ( i=0; i<100; i++ ) { + + signed int h, k, l; + double x, y, z; + double nh, nk, nl; + double dh, dk, dl; + int f = 0; + long int ih, ik, il; + + h = flat_noise(0, 30); + k = flat_noise(0, 30); + l = flat_noise(0, 30); + + x = h*asx + k*bsx + l*csx; + y = h*asy + k*bsy + l*csy; + z = h*asz + k*bsz + l*csz; + + nh = x*ax + y*ay + z*az; + nk = x*bx + y*by + z*bz; + nl = x*cx + y*cy + z*cz; + + dh = nh - lrint(nh); dk = nk - lrint(nk); dl = nl - lrint(nl); + + if ( fabs(dh) > 0.1 ) f++; + if ( fabs(dk) > 0.1 ) f++; + if ( fabs(dl) > 0.1 ) f++; + + ih = lrint(nh); ik = lrint(nk); il = lrint(nl); + if ( forbidden_reflection(cell, ih, ik, il) ) { + STATUS("Primitive %3i %3i %3i -> " + "Centered %3li %3li %3li, " + "which is forbidden\n", + h, k, l, ih, ik, il); + fail = 1; + } + + if ( f ) { + STATUS("Primitive %3i %3i %3i -> " + "Centered %7.2f %7.2f %7.2f\n", + h, k, l, nh, nk, nl); + fail = 1; + } + + } return fail; } - int main(int argc, char *argv[]) { int fail = 0; @@ -133,15 +230,15 @@ int main(int argc, char *argv[]) /* Orthorhombic A */ fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, - L_ORTHORHOMBIC, 'A', '*'); + L_ORTHORHOMBIC, 'A', 'a'); /* Orthorhombic B */ fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, - L_ORTHORHOMBIC, 'B', '*'); + L_ORTHORHOMBIC, 'B', 'b'); /* Orthorhombic C */ fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, - L_ORTHORHOMBIC, 'C', '*'); + L_ORTHORHOMBIC, 'C', 'c'); /* Orthorhombic I */ fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0, @@ -180,10 +277,6 @@ int main(int argc, char *argv[]) L_HEXAGONAL, 'P', 'c'); /* Hexagonal H (PDB-speak for rhombohedral) */ - fail += check_centering(40e-10, 20e-10, 20e-10, 120.0, 90.0, 90.0, - L_HEXAGONAL, 'H', 'a'); - fail += check_centering(20e-10, 40e-10, 20e-10, 90.0, 120.0, 90.0, - L_HEXAGONAL, 'H', 'b'); fail += check_centering(20e-10, 20e-10, 40e-10, 90.0, 90.0, 120.0, L_HEXAGONAL, 'H', 'c'); diff --git a/tests/transformation_check.c b/tests/transformation_check.c index 80324795..7d25aa04 100644 --- a/tests/transformation_check.c +++ b/tests/transformation_check.c @@ -73,6 +73,39 @@ static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn) } +static int check_identity(UnitCell *cell, UnitCellTransformation *tfn) +{ + UnitCell *cnew; + double a[9], b[9]; + int i; + int fail = 0; + + cnew = cell_transform(cell, tfn); + + cell_get_cartesian(cell, &a[0], &a[1], &a[2], + &a[3], &a[4], &a[5], + &a[6], &a[7], &a[8]); + cell_get_cartesian(cnew, &b[0], &b[1], &b[2], + &b[3], &b[4], &b[5], + &b[6], &b[7], &b[8]); + for ( i=0; i<9; i++ ) { + if ( !within_tolerance(a[i], b[i], 0.1) ) { + fail = 1; + STATUS("%e %e\n", a[i], b[i]); + } + } + + if ( fail ) { + ERROR("Original cell not recovered after transformation:\n"); + cell_print(cell); + tfn_print(tfn); + cell_print(cnew); + } + + return fail; +} + + int main(int argc, char *argv[]) { int fail = 0; @@ -125,6 +158,18 @@ int main(int argc, char *argv[]) fail += check_transformation(cell, tfn); tfn_free(tfn); + /* Identity in two parts */ + tfn = tfn_identity(); + if ( tfn == NULL ) return 1; + tfn_combine(tfn, tfn_vector(0,0,1), + tfn_vector(0,1,0), + tfn_vector(-1,0,0)); + tfn_combine(tfn, tfn_vector(0,0,-1), + tfn_vector(0,1,0), + tfn_vector(1,0,0)); + fail += check_identity(cell, tfn); + tfn_free(tfn); + cell_free(cell); return fail; -- cgit v1.2.3 From 9aafac1398cfe51c3082f3a708c1fe044f0f6884 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 1 Oct 2012 18:37:48 +0200 Subject: Remove debugging message --- libcrystfel/src/cell-utils.c | 3 --- 1 file changed, 3 deletions(-) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index f854efe0..3aaae90e 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -989,9 +989,6 @@ UnitCell *load_cell_from_pdb(const char *filename) 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); -- cgit v1.2.3 From d51007fe1a9d8a797bdd9b71eea715a09252cc6e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 1 Oct 2012 17:05:43 +0200 Subject: render_hkl: New option "--zone" for selecting Laue zone number --- src/render_hkl.c | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/render_hkl.c b/src/render_hkl.c index 3ae3b345..7818b304 100644 --- a/src/render_hkl.c +++ b/src/render_hkl.c @@ -158,7 +158,7 @@ static void draw_circles(double xh, double xk, double xl, cairo_t *dctx, int wght, double boost, int colscale, UnitCell *cell, double radius, double theta, double as, double bs, double cx, double cy, - double scale, double max_val) + double scale, double max_val, signed int zone) { Reflection *refl; RefListIterator *iter; @@ -193,7 +193,7 @@ static void draw_circles(double xh, double xk, double xl, get_equiv(sym, m, i, ha, ka, la, &h, &k, &l); /* Is the reflection in the zone? */ - if ( h*zh + k*zk + l*zl != 0 ) continue; + if ( h*zh + k*zk + l*zl != zone ) continue; xi = (h*xh + k*xk + l*xl) / nx; yi = (h*yh + k*yk + l*yl) / ny; @@ -300,7 +300,7 @@ static void render_za(UnitCell *cell, RefList *list, int colscale, signed int xh, signed int xk, signed int xl, signed int yh, signed int yk, signed int yl, - const char *outfile, double scale_top) + const char *outfile, double scale_top, signed int zone) { cairo_surface_t *surface; cairo_t *dctx; @@ -416,7 +416,7 @@ static void render_za(UnitCell *cell, RefList *list, draw_circles(xh, xk, xl, yh, yk, yl, zh, zk, zl, list, sym, dctx, wght, boost, colscale, cell, max_r, theta, as, bs, cx, cy, scale, - max_val); + max_val, zone); /* Centre marker */ cairo_arc(dctx, (double)cx, @@ -610,6 +610,7 @@ int main(int argc, char *argv[]) char *outfile = NULL; double scale_top = -1.0; char *endptr; + long int zone = 0; /* Long options */ const struct option longopts[] = { @@ -626,6 +627,7 @@ int main(int argc, char *argv[]) {"counts", 0, &config_sqrt, 1}, {"colour-key", 0, &config_colkey, 1}, {"scale-top", 1, NULL, 2}, + {"zone", 1, NULL, 3}, {0, 0, NULL, 0} }; @@ -686,6 +688,17 @@ int main(int argc, char *argv[]) } break; + case 3 : + errno = 0; + zone = strtol(optarg, &endptr, 10); + if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) + || (errno != 0) ) + { + ERROR("Invalid zone number ('%s')\n", optarg); + return 1; + } + break; + case 0 : break; @@ -805,7 +818,7 @@ int main(int argc, char *argv[]) } render_za(cell, list, boost, sym, wght, colscale, - rh, rk, rl, dh, dk, dl, outfile, scale_top); + rh, rk, rl, dh, dk, dl, outfile, scale_top, zone); free(pdb); free_symoplist(sym); -- cgit v1.2.3