Finishing tweaks for uncentering stuff

Today, I have mostly been having my life made difficult by the PDB's invention of "H centering".

2 files changed, 171 insertions, 33 deletions

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;