Tighter symmetry check

src/check_hkl.c

@@ -379,6 +379,13 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	/* Check that the intensities have the correct symmetry */
+	if ( check_symmetry(items, sym) ) {
+		ERROR("The input reflection list does not appear to"
+		      " have symmetry %s\n", sym);
+		return 1;
+	}
+
 	/* Reject reflections */
 	good_items = new_items();
 	for ( i=0; i<num_items(items); i++ ) {

src/compare_hkl.c

@@ -352,6 +352,18 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	/* Check that the intensities have the correct symmetry */
+	if ( check_symmetry(i1, sym) ) {
+		ERROR("The first input reflection list does not appear to"
+		      " have symmetry %s\n", sym);
+		return 1;
+	}
+	if ( check_symmetry(i2, sym) ) {
+		ERROR("The second input reflection list does not appear to"
+		      " have symmetry %s\n", sym);
+		return 1;
+	}
+
 	/* List for output scale factor map */
 	out = new_list_intensity();
 

src/get_hkl.c

@@ -358,6 +358,11 @@ int main(int argc, char *argv[])
 		input_items = read_reflections(input, ideal_ref, phases,
 		                               NULL, esds);
 		free(input);
+		if ( check_symmetry(input_items, mero) ) {
+			ERROR("The input reflection list does not appear to"
+			      " have symmetry %s\n", mero);
+			return 1;
+		}
 	}
 
 	if ( config_poisson ) poisson_reflections(ideal_ref, input_items);

src/indexamajig.c

@@ -745,6 +745,12 @@ int main(int argc, char *argv[])
 			set_flag(flags, it->h, it->k, it->l, 1);
 		}
 
+		if ( check_symmetry(items, sym) ) {
+			ERROR("The input reflection list does not appear to"
+			      " have symmetry %s\n", sym);
+			return 1;
+		}
+
 		delete_items(items);
 
 	} else {

src/pattern_sim.c

@@ -32,6 +32,7 @@
 #include "sfac.h"
 #include "reflections.h"
 #include "beam-parameters.h"
+#include "symmetry.h"
 
 
 static void show_help(const char *s)
@@ -445,6 +446,13 @@ int main(int argc, char *argv[])
 			set_flag(flags, it->h, it->k, it->l, 1);
 		}
 
+		/* Check that the intensities have the correct symmetry */
+		if ( check_symmetry(items, sym) ) {
+			ERROR("The input reflection list does not appear to"
+			      " have symmetry %s\n", sym);
+			return 1;
+		}
+
 		delete_items(items);
 	}
 

src/render_hkl.c

@@ -648,6 +648,11 @@ int main(int argc, char *argv[])
 		ERROR("Couldn't open file '%s'\n", infile);
 		return 1;
 	}
+	if ( check_symmetry(items, sym) ) {
+		ERROR("The input reflection list does not appear to"
+		      " have symmetry %s\n", sym);
+		return 1;
+	}
 
 	if ( config_povray ) {
 		r = povray_render_animation(cell, ref, cts, items,

src/symmetry.c

@@ -454,22 +454,6 @@ ReflItemList *get_twins(ReflItemList *items, const char *holo, const char *mero)
 }
 
 
-static void scold_user_about_symmetry(signed int h, signed int k, signed int l,
-                                      signed int he, signed int ke,
-                                      signed int le)
-{
-	ERROR("Symmetrically equivalent reflection (%i %i %i) found for "
-	      "%i %i %i in the input.\n", he, ke, le, h, k, l);
-	ERROR("This indicates that you lied to me about the symmetry of the "
-	      "input reflections.  ");
-	ERROR("I won't be able to give you a meaningful result in this "
-	      "situation, so I'm going to give up right now.  ");
-	ERROR("Please reconsider your previous processing of the data, and "
-	      "perhaps try again with a lower symmetry for the '-y' option.\n");
-	abort();
-}
-
-
 int find_unique_equiv(ReflItemList *items, signed int h, signed int k,
                       signed int l, const char *mero, signed int *hu,
                       signed int *ku, signed int *lu)
@@ -485,19 +469,18 @@ int find_unique_equiv(ReflItemList *items, signed int h, signed int k,
 		f = find_item(items, he, ke, le);
 
 		/* There must only be one equivalent.  If there are more, it
-		 * indicates that the user lied about the input symmetry. */
-		if ( f && found ) {
-			scold_user_about_symmetry(he, ke, le, *hu, *ku, *lu);
-		}
+		 * indicates that the user lied about the input symmetry.
+		 * This situation should have been checked for earlier by
+		 * calling check_symmetry() with 'items' and 'mero'. */
 
 		if ( f && !found ) {
 			*hu = he;  *ku = ke;  *lu = le;
-			found = 1;
+			return 1;
 		}
 
 	}
 
-	return found;
+	return 0;
 }
 
 
@@ -615,3 +598,44 @@ int has_bisecting_mirror_or_diad(const char *sym)
 	ERROR("Couldn't find mirror definition for '%s'.\n", sym);
 	abort();
 }
+
+
+int check_symmetry(ReflItemList *items, const char *sym)
+{
+	int i;
+	unsigned char *flags;
+
+	flags = new_list_flag();
+	for ( i=0; i<num_items(items); i++ ) {
+		struct refl_item *it = get_item(items, i);
+		set_flag(flags, it->h, it->k, it->l, 1);
+	}
+
+	for ( i=0; i<num_items(items); i++ ) {
+
+		int j;
+		struct refl_item *it = get_item(items, i);
+		int found = 0;
+
+		for ( j=0; j<num_equivs(it->h, it->k, it->l, sym); j++ ) {
+
+			signed int he, ke, le;
+			get_equiv(it->h, it->k, it->l, &he, &ke, &le, sym, j);
+
+			if ( abs(he) > INDMAX ) continue;
+			if ( abs(le) > INDMAX ) continue;
+			if ( abs(ke) > INDMAX ) continue;
+
+			found += lookup_flag(flags, he, ke, le);
+
+		}
+
+		if ( found > 1 ) {
+			STATUS("%i found for %i %i %i\n", found, it->h, it->k, it->l);
+			return 1;  /* Symmetry is wrong! */
+		}
+
+	}
+
+	return 0;
+}

src/symmetry.h

@@ -42,6 +42,8 @@ extern int find_unique_equiv(ReflItemList *items, signed int h, signed int k,
                              signed int l, const char *mero, signed int *hu,
                              signed int *ku, signed int *lu);
 
+extern int check_symmetry(ReflItemList *items, const char *sym);
+
 /* Properties of point groups */
 extern int is_polyhedral(const char *sym);
 extern int rotational_order(const char *sym);