From 5d0cae2baa9a7e09dd54b144cc503feac730794a Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Tue, 13 Jul 2010 11:57:21 +0200
Subject: Simplify symmetry and twinning quite a lot

---
 src/symmetry.c | 229 +++++++++++++++++++++++++++++++--------------------------
 src/symmetry.h |  10 ++-
 src/utils.c    |  21 +++++-
 src/utils.h    |   3 +
 4 files changed, 152 insertions(+), 111 deletions(-)

(limited to 'src')

diff --git a/src/symmetry.c b/src/symmetry.c
index e0b92335..b29af3d4 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -57,35 +57,16 @@ static int check_cond(signed int h, signed int k, signed int l, const char *sym)
 	}
 
 
-int num_equivs(signed int h, signed int k, signed int l, const char *sym)
+static int num_general_equivs(const char *sym)
 {
-	/* 000 is always unique */
-	if ( (h==0) && (k==0) && (l==0) ) return 1;
-
+	/* Triclinic */
 	if ( strcmp(sym, "1") == 0 ) return 1;
-
 	if ( strcmp(sym, "-1") == 0 ) return 2;
 
-	if ( strcmp(sym, "6") == 0 ) {
-		if ( (h==0) && (k==0) ) return 1;  /* a */
-		return 6;  /* b */
-	}
-
-	if ( strcmp(sym, "6/m") == 0 ) {
-		if ( (h==0) && (k==0) ) return 2;  /* a */
-		if ( l == 0 ) return 6;  /* b */
-		return 12;  /* c */
-	}
-
-	if ( strcmp(sym, "6/mmm") == 0 ) {
-		if ( (h==0) && (k==0) ) return 2;  /* a */
-		if ( (l==0) && (h==k) ) return 6;  /* b */
-		if ( (l==0) && (k==0) ) return 6;  /* c */
-		if ( h == k ) return 12;  /* d */
-		if ( k == 0 ) return 12;  /* e */
-		if ( l == 0 ) return 12;  /* f */
-		return 24;  /* g */
-	}
+	/* Hexagonal */
+	if ( strcmp(sym, "6") == 0 ) return 6;
+	if ( strcmp(sym, "6/m") == 0 ) return 12;
+	if ( strcmp(sym, "6/mmm") == 0 ) return 24;
 
 	/* TODO: Add more groups here */
 
@@ -93,9 +74,9 @@ int num_equivs(signed int h, signed int k, signed int l, const char *sym)
 }
 
 
-void get_equiv(signed int h, signed int k, signed int l,
-               signed int *he, signed int *ke, signed int *le,
-               const char *sym, int idx)
+void get_general_equiv(signed int h, signed int k, signed int l,
+                       signed int *he, signed int *ke, signed int *le,
+                       const char *sym, int idx)
 {
 	signed int i = -h-k;
 
@@ -111,13 +92,6 @@ void get_equiv(signed int h, signed int k, signed int l,
 	}
 
 	if ( strcmp(sym, "6") == 0 ) {
-		/* a */
-		if ( (h==0) && (k==0) ) {
-			switch ( idx ) {
-			case  0 : *he = 0;   *ke = 0;   *le = l;   return;
-			}
-		}
-		/* b */
 		switch ( idx ) {
 		case  0 : *he = h;   *ke = k;   *le = l;  return;
 		case  1 : *he = i;   *ke = h;   *le = l;  return;
@@ -129,14 +103,6 @@ void get_equiv(signed int h, signed int k, signed int l,
 	}
 
 	if ( strcmp(sym, "6/m") == 0 ) {
-		/* a */
-		if ( (h==0) && (k==0) ) {
-			switch ( idx ) {
-			case  0 : *he = 0;   *ke = 0;   *le = l;   return;
-			case  1 : *he = 0;   *ke = 0;   *le = -l;  return;
-			}
-		}
-		/* b-c */
 		switch ( idx ) {
 		case  0 : *he = h;   *ke = k;   *le = l;   return;
 		case  1 : *he = i;   *ke = h;   *le = l;   return;
@@ -154,14 +120,6 @@ void get_equiv(signed int h, signed int k, signed int l,
 	}
 
 	if ( strcmp(sym, "6/mmm") == 0 ) {
-		/* a */
-		if ( (h==0) && (k==0) ) {
-			switch ( idx ) {
-			case  0 : *he = 0;   *ke = 0;   *le = l;   return;
-			case  1 : *he = 0;   *ke = 0;   *le = -l;  return;
-			}
-		}
-		/* b-g */
 		switch ( idx ) {
 		case  0 : *he = h;   *ke = k;   *le = l;   return;
 		case  1 : *he = i;   *ke = h;   *le = l;   return;
@@ -197,6 +155,66 @@ void get_equiv(signed int h, signed int k, signed int l,
 }
 
 
+static int special_position(signed int hs, signed int ks, signed int ls,
+                            signed int *hp, signed int *kp, signed int *lp,
+                            const char *sym, signed int idx)
+{
+	int n_general;
+	int i;
+	ReflItemList *equivs;
+	int n_equivs = 0;
+
+	equivs = new_items();
+
+	if ( idx == 0 ) {
+		/* Index zero is always the original reflection */
+		*hp = hs;  *kp = ks;  *lp = ls;
+		return 0;
+	}
+
+	n_general = num_general_equivs(sym);
+
+	for ( i=0; i<n_general; i++ ) {
+
+		signed int h, k, l;
+
+		/* Get equivalent according to the holohedral group */
+		get_general_equiv(hs, ks, ls, &h, &k, &l, sym, i);
+
+		/* Already got this one? */
+		if ( find_item(equivs, h, k, l) ) continue;
+
+		if ( n_equivs == idx ) {
+			*hp = h;
+			*kp = k;
+			*lp = l;
+			delete_items(equivs);
+			return n_equivs;
+		}
+		add_item(equivs, h, k, l);
+		n_equivs++;
+
+	}
+
+	delete_items(equivs);
+	return n_equivs;
+}
+
+
+void get_equiv(signed int h, signed int k, signed int l,
+               signed int *he, signed int *ke, signed int *le,
+               const char *sym, int idx)
+{
+	special_position(h, k, l, he, ke, le, sym, idx);
+}
+
+
+int num_equivs(signed int h, signed int k, signed int l, const char *sym)
+{
+	return special_position(h, k, l, NULL, NULL, NULL, sym, -1);
+}
+
+
 void get_asymm(signed int h, signed int k, signed int l,
                signed int *hp, signed int *kp, signed int *lp,
                const char *sym)
@@ -218,7 +236,7 @@ void get_asymm(signed int h, signed int k, signed int l,
 }
 
 
-static const char *get_holohedral(const char *sym)
+const char *get_holohedral(const char *sym)
 {
 	/* Triclinic */
 	if ( strcmp(sym, "1") == 0 ) return "-1";
@@ -237,37 +255,28 @@ static const char *get_holohedral(const char *sym)
 
 
 /* This is kind of like a "numerical" left coset decomposition.
- * Given a reflection index and a point group, it returns the "idx-th"
- * twinning possibility for the reflection.  It returns "idx" if
- * successful.  To just count the number of possibilities, set idx=-1.
+ * Given a reflection index and a point group, it returns the list of twinning
+ * possibilities.
  *
- * The sequence of operators producing each possibility is guaranteed to
- * be the same for any choice of indices given the same point group. */
-static int coset_decomp(signed int hs, signed int ks, signed int ls,
-                        signed int *hp, signed int *kp, signed int *lp,
-                        const char *mero, signed int idx)
+ * To count the number of possibilities, use num_items() on the result.
+ */
+static ReflItemList *coset_decomp(signed int hs, signed int ks, signed int ls,
+                                  const char *mero)
 {
 	const char *holo = get_holohedral(mero);
 	int n_mero, n_holo;
 	int i;
-	signed int n_twins = 1;
 	signed int h, k, l;
-	ReflItemList *twins;
-
-	twins = new_items();
-
-	if ( idx == 0 ) {
-		/* Twin index zero is always the original orientation */
-		*hp = hs;  *kp = ks;  *lp = ls;
-		return 0;
-	}
+	ReflItemList *twins = new_items();
 
+	/* Start by putting the given reflection into the asymmetric cell
+	 * for its (probably merohedral) point group. */
 	get_asymm(hs, ks, ls, &h, &k, &l, mero);
 
 	/* How many equivalents in the holohedral point group are not
 	 * equivalents according to the (possibly) merohedral group? */
-	n_holo = num_equivs(h, k, l, holo);
-	n_mero = num_equivs(h, k, l, mero);
+	n_holo = num_general_equivs(holo);
+	n_mero = num_general_equivs(mero);
 
 	for ( i=0; i<n_holo; i++ ) {
 
@@ -275,51 +284,61 @@ static int coset_decomp(signed int hs, signed int ks, signed int ls,
 		signed int hs_holo, ks_holo, ls_holo;
 
 		/* Get equivalent according to the holohedral group */
-		get_equiv(h, k, l, &hs_holo, &ks_holo, &ls_holo, holo, i);
+		get_general_equiv(h, k, l, &hs_holo, &ks_holo, &ls_holo,
+		                  holo, i);
 
 		/* Put it into the asymmetric cell for the merohedral group */
 		get_asymm(hs_holo, ks_holo, ls_holo,
 		          &h_holo, &k_holo, &l_holo, mero);
 
-		/* Is this the same reflection as we started with?
-		 * If not, this reflection is 'equivalent by twinning' */
-		if ( (h_holo != h) || (k_holo != k) || (l_holo != l) ) {
-
-			if ( find_item(twins, h_holo, k_holo, l_holo) )
-				continue;
-
-			if ( n_twins == idx ) {
-				*hp = h_holo;
-				*kp = k_holo;
-				*lp = l_holo;
-				delete_items(twins);
-				return n_twins;
-			}
-			add_item(twins, h_holo, k_holo, l_holo);
-			n_twins++;
-		}
+		/* Already got this one?
+		 * Note: The list "twins" starts empty, so the first iteration
+		 * (i=0) will add the original reflection to the list along with
+		 * the identity operation. */
+		if ( find_item(twins, h_holo, k_holo, l_holo) ) continue;
+
+		add_item_with_op(twins, h_holo, k_holo, l_holo, i);
 
 	}
 
-	delete_items(twins);
-	return n_twins;
+	return twins;
 }
 
 
-/* Get the number of twinned "equivalents" for this reflection */
-int num_twins(signed int h, signed int k, signed int l, const char *sym)
+/* Work out the twinning possibilities for this pattern.
+ * To use the result, call get_general_equiv() on each reflection using
+ * the holohedral point group (use get_holohedral() for this), and for "idx"
+ * give each "op" field from the list returned by this function. */
+ReflItemList *get_twins(ReflItemList *items, const char *sym)
 {
-	return coset_decomp(h, k, l, NULL, NULL, NULL, sym, -1);
-}
-
+	int i;
+	int n_twins = 1;
+	ReflItemList *max_ops = NULL;
+
+	/* Run the coset decomposition for every reflection in the "pattern",
+	 * and see which gives the highest number of possibilities.  This
+	 * approach takes into account that a pattern consisting entirely of
+	 * special reflections might have fewer twin possibilities. */
+	for ( i=0; i<num_items(items); i++ ) {
+
+		signed int h, k, l;
+		struct refl_item *item;
+		ReflItemList *ops;
+
+		item = get_item(items, i);
+
+		h = item->h;
+		k = item->k;
+		l = item->l;
+
+		ops = coset_decomp(h, k, l, sym);
+		if ( num_items(ops) > n_twins ) {
+			n_twins = num_items(ops);
+			delete_items(max_ops);
+			max_ops = ops;
+		}
 
-void get_twins(signed int h, signed int k, signed int l,
-              signed int *hp, signed int *kp, signed int *lp,
-              const char *sym, int idx)
-{
-	if ( coset_decomp(h, k, l, hp, kp, lp, sym, idx) != idx ) {
-		ERROR("Failed coset decomposition for %i %i %i in %s\n",
-		      h, k, l, sym);
-		abort();
 	}
+
+	return max_ops;
 }
diff --git a/src/symmetry.h b/src/symmetry.h
index 52dd387d..0ae01e53 100644
--- a/src/symmetry.h
+++ b/src/symmetry.h
@@ -29,11 +29,13 @@ extern void get_equiv(signed int h, signed int k, signed int l,
                       signed int *he, signed int *ke, signed int *le,
                       const char *sym, int idx);
 
-extern int num_twins(signed int h, signed int k, signed int l, const char *sym);
+extern void get_general_equiv(signed int h, signed int k, signed int l,
+                              signed int *he, signed int *ke, signed int *le,
+                              const char *sym, int idx);
 
-extern void get_twins(signed int h, signed int k, signed int l,
-                      signed int *hp, signed int *kp, signed int *lp,
-                      const char *sym, int idx);
+extern const char *get_holohedral(const char *sym);
+
+extern ReflItemList *get_twins(ReflItemList *items, const char *sym);
 
 
 #endif	/* SYMMETRY_H */
diff --git a/src/utils.c b/src/utils.c
index 649d2b1c..1e54f923 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -292,17 +292,20 @@ struct _reflitemlist {
 	int max_items;
 };
 
+
 void clear_items(ReflItemList *items)
 {
 	items->n_items = 0;
 }
 
+
 static void alloc_items(ReflItemList *items)
 {
 	items->items = realloc(items->items,
 	                       items->max_items*sizeof(struct refl_item));
 }
 
+
 ReflItemList *new_items()
 {
 	ReflItemList *new;
@@ -314,14 +317,17 @@ ReflItemList *new_items()
 	return new;
 }
 
+
 void delete_items(ReflItemList *items)
 {
+	if ( items == NULL ) return;
 	if ( items->items != NULL ) free(items->items);
 	free(items);
 }
 
-void add_item(ReflItemList *items,
-                     signed int h, signed int k, signed int l)
+
+void add_item_with_op(ReflItemList *items, signed int h, signed int k,
+                      signed int l, int op)
 {
 	if ( items->n_items == items->max_items ) {
 		items->max_items += 1024;
@@ -331,9 +337,17 @@ void add_item(ReflItemList *items,
 	items->items[items->n_items].h = h;
 	items->items[items->n_items].k = k;
 	items->items[items->n_items].l = l;
+	items->items[items->n_items].op = op;
 	items->n_items++;
 }
 
+
+void add_item(ReflItemList *items, signed int h, signed int k, signed int l)
+{
+	add_item_with_op(items, h, k, l, 0);
+}
+
+
 int find_item(ReflItemList *items,
                      signed int h, signed int k, signed int l)
 {
@@ -348,17 +362,20 @@ int find_item(ReflItemList *items,
 	return 0;
 }
 
+
 struct refl_item *get_item(ReflItemList *items, int i)
 {
 	if ( i >= items->n_items ) return NULL;
 	return &items->items[i];
 }
 
+
 int num_items(const ReflItemList *items)
 {
 	return items->n_items;
 }
 
+
 unsigned int *items_to_counts(ReflItemList *items)
 {
 	unsigned int *c;
diff --git a/src/utils.h b/src/utils.h
index ab07336d..46ce51a6 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -184,6 +184,7 @@ struct refl_item {
 	signed int h;
 	signed int k;
 	signed int l;
+	int op;
 };
 
 extern void clear_items(ReflItemList *items);
@@ -191,6 +192,8 @@ extern ReflItemList *new_items(void);
 extern void delete_items(ReflItemList *items);
 extern void add_item(ReflItemList *items,
                      signed int h, signed int k, signed int l);
+extern void add_item_with_op(ReflItemList *items,
+                             signed int h, signed int k, signed int l, int op);
 extern int find_item(ReflItemList *items,
                      signed int h, signed int k, signed int l);
 extern struct refl_item *get_item(ReflItemList *items, int i);
-- 
cgit v1.2.3