diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/symmetry.c | 229 | ||||
-rw-r--r-- | src/symmetry.h | 10 | ||||
-rw-r--r-- | src/utils.c | 21 | ||||
-rw-r--r-- | src/utils.h | 3 |
4 files changed, 152 insertions, 111 deletions
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); |