aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/cell-utils.c128
-rw-r--r--libcrystfel/src/cell-utils.h5
-rw-r--r--libcrystfel/src/cell.c11
-rw-r--r--tests/centering_check.c159
-rw-r--r--tests/transformation_check.c45
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 <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
+#include <fenv.h>
#include <cell.h>
#include <cell-utils.h>
@@ -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;
generated by cgit v1.2.3 (git 2.25.1) at 2025-01-22 02:43:55 +0000