diff options
| -rw-r--r-- | src/process_hkl.c | 84 | ||||
| -rw-r--r-- | src/symmetry.c | 7 | ||||
| -rw-r--r-- | src/symmetry.h | 3 |
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 */ |