process_hkl: Improve de-twinning

3 files changed, 71 insertions, 23 deletions

diff --git a/src/process_hkl.c b/src/process_hkl.c
index 3a9d8d34..29288154 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -180,21 +180,49 @@ static void process_reflections(double *ref, unsigned int *counts,
 }
 
 
-static int resolve_twin(ReflItemList *items, const char *symm,
-                        const double *patt, const double *model,
-                        const unsigned int *model_counts)
+/* Note "holo" needn't actually be a holohedral point group, if you want to try
+ * something strange like resolving from a low-symmetry group into an even
+ * lower symmetry one.
+ */
+static ReflItemList *get_twin_possibilities(const char *holo, const char *mero)
 {
+	ReflItemList *test_items;
 	ReflItemList *twins;
-	const char *holo;
+	int np;
+
+	np = num_general_equivs(holo) / num_general_equivs(mero);
+
+	test_items = new_items();
+
+	/* Some arbitrarily chosen reflections which can't be special
+	 * reflections in any point group, i.e. lots of odd numbers,
+	 * prime numbers and so on.  There's probably an analytical
+	 * way of working these out, but this will do. */
+	add_item(test_items, 1, 2, 3);
+	add_item(test_items, 3, 7, 13);
+	add_item(test_items, 5, 2, 1);
+
+	twins = get_twins(test_items, holo, mero);
+	delete_items(test_items);
+
+	if ( num_items(twins) != np ) {
+		ERROR("Whoops! Couldn't find all the twinning possiblities.\n");
+		abort();
+	}
+
+	return twins;
+}
+
+
+static int resolve_twin(ReflItemList *items, ReflItemList *twins,
+                        const double *patt, const double *model,
+                        const unsigned int *model_counts,
+                        const char *holo, const char *mero)
+{
 	int n, i;
 	double best_fom = 0.0;
 	int best_op = 0;
 
-	/* How many possible twinned orientations are there? */
-	twins = get_twins(items, symm);
-
-	holo = get_holohedral(symm);
-
 	n = num_items(twins);
 
 	for ( i=0; i<n; i++ ) {
@@ -214,7 +242,7 @@ static int resolve_twin(ReflItemList *items, const char *symm,
 
 			get_general_equiv(r->h, r->k, r->l, &h, &k, &l,
 			                  holo, op);
-			get_asymm(h, k, l, &h, &k, &l, symm);
+			get_asymm(h, k, l, &h, &k, &l, mero);
 
 			set_intensity(trial_ints, h, k, l,
 			              lookup_intensity(patt, r->h, r->k, r->l));
@@ -237,7 +265,6 @@ static int resolve_twin(ReflItemList *items, const char *symm,
 	}
 	//printf("\n");
 
-	delete_items(twins);
 	return best_op;
 }
 
@@ -245,13 +272,18 @@ static int resolve_twin(ReflItemList *items, const char *symm,
 static void merge_pattern(double *model, const double *new,
                           unsigned int *model_counts,
                           ReflItemList *items, int mo, int sum,
-                          const char *symm)
+                          ReflItemList *twins,
+                          const char *holo, const char *mero)
 {
 	int i;
-	const char *holo = get_holohedral(symm);  /* Needed to look up later */
 	int twin;
 
-	twin = resolve_twin(items, symm, new, model, model_counts);
+	if ( twins != NULL ) {
+		twin = resolve_twin(items, twins, new, model,
+		                     model_counts, holo, mero);
+	} else {
+		twin = 0;
+	}
 
 	for ( i=0; i<num_items(items); i++ ) {
 
@@ -267,7 +299,7 @@ static void merge_pattern(double *model, const double *new,
 		ls = item->l;
 
 		get_general_equiv(hs, ks, ls, &h, &k, &l, holo, twin);
-		get_asymm(h, k, l, &h, &k, &l, symm);
+		get_asymm(h, k, l, &h, &k, &l, mero);
 
 		intensity = lookup_intensity(new, h, k, l);
 
@@ -317,6 +349,8 @@ int main(int argc, char *argv[])
 	int f0_valid;
 	int n_nof0 = 0;
 	ReflItemList *items;
+	ReflItemList *twins;
+	int i;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -445,9 +479,22 @@ int main(int argc, char *argv[])
 	/* Show useful symmetry information */
 	const char *holo = get_holohedral(sym);
 	int np = num_general_equivs(holo) / num_general_equivs(sym);
-	STATUS("Resolving from point group %s to %s (%i possibilities)\n",
-	       holo, sym, np);
+	if ( np > 1 ) {
+
+		STATUS("Resolving point group %s into %s (%i possibilities)\n",
+		       holo, sym, np);
+		/* Get the list of twin/Bijvoet possibilities */
+		twins = get_twin_possibilities(holo, sym);
+		STATUS("Twin/inversion operation indices (from %s) are:", holo);
+		for ( i=0; i<num_items(twins); i++ ) {
+			STATUS(" %i", get_item(twins, i)->op);
+		}
+		STATUS("\n");
 
+	} else {
+		STATUS("No twin/inversion resolution necessary.\n");
+		twins = NULL;
+	}
 
 	n_patterns = 0;
 	f0_valid = 0;
@@ -489,7 +536,8 @@ int main(int argc, char *argv[])
 
 			/* Start of second or later pattern */
 			merge_pattern(model, new_pattern, model_counts,
-			              items, config_maxonly, config_sum, sym);
+			              items, config_maxonly, config_sum, twins,
+			              holo, sym);
 
 			if ( (trueref != NULL) && config_every
 			    && (n_patterns % config_every == 0) ) {
diff --git a/src/symmetry.c b/src/symmetry.c
index 31470db7..69a2e55c 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -265,9 +265,8 @@ const char *get_holohedral(const char *sym)
  * 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, const char *mero)
 {
-	const char *holo = get_holohedral(mero);
 	int n_mero, n_holo;
 	int i;
 	signed int h, k, l;
@@ -313,7 +312,7 @@ static ReflItemList *coset_decomp(signed int hs, signed int ks, signed int ls,
  * 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)
+ReflItemList *get_twins(ReflItemList *items, const char *holo, const char *mero)
 {
 	int i;
 	ReflItemList *ops = new_items();;
@@ -334,7 +333,7 @@ ReflItemList *get_twins(ReflItemList *items, const char *sym)
 		k = item->k;
 		l = item->l;
 
-		new_ops = coset_decomp(h, k, l, sym);
+		new_ops = coset_decomp(h, k, l, holo, mero);
 		union_op_items(ops, new_ops);
 		delete_items(new_ops);
 
diff --git a/src/symmetry.h b/src/symmetry.h
index 6173166c..437066ad 100644
--- a/src/symmetry.h
+++ b/src/symmetry.h
@@ -37,7 +37,8 @@ extern void get_general_equiv(signed int h, signed int k, signed int l,
 
 extern const char *get_holohedral(const char *sym);
 
-extern ReflItemList *get_twins(ReflItemList *items, const char *sym);
+extern ReflItemList *get_twins(ReflItemList *items,
+                               const char *holo, const char *mero);
 
 
 #endif	/* SYMMETRY_H */