aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--doc/reference/CrystFEL-sections.txt17
-rw-r--r--src/symmetry.c868
-rw-r--r--src/symmetry.h14
3 files changed, 274 insertions, 625 deletions
diff --git a/doc/reference/CrystFEL-sections.txt b/doc/reference/CrystFEL-sections.txt
index c032fe76..b9135f1b 100644
--- a/doc/reference/CrystFEL-sections.txt
+++ b/doc/reference/CrystFEL-sections.txt
@@ -195,17 +195,14 @@ indexingprivate
<SECTION>
<FILE>symmetry</FILE>
-check_list_symmetry
-check_symmetry
-has_bisecting_mirror_or_diad
-has_perpendicular_mirror
-is_polyhedral
-rotational_order
-num_equivs
-num_general_equivs
-get_asymm
+SymOpList
+<SUBSECTION>
+free_symoplist
+num_ops
+get_pointgroup
get_equiv
-get_general_equiv
+special_position
+get_asymm
get_twins
</SECTION>
diff --git a/src/symmetry.c b/src/symmetry.c
index 3a9d505e..77ac3ec3 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -23,12 +23,6 @@
#include "utils.h"
-#ifdef DEBUG
-#define SYM_DEBUG STATUS
-#else /* DEBUG */
-#define SYM_DEBUG(...)
-#endif /* DEBUG */
-
/**
* SECTION:symmetry
* @short_description: Point symmetry handling
@@ -42,11 +36,26 @@
*/
-struct sym_op {
- signed int h;
- signed int k;
- signed int l;
- int op;
+enum lattice_type
+{
+ L_TRICLINIC,
+ L_MONOCLINIC,
+ L_ORTHORHOMBIC,
+ L_TETRAGONAL,
+ L_RHOMBOHEDRAL,
+ L_TRIGONAL,
+ L_HEXAGONAL,
+ L_CUBIC,
+};
+
+
+struct sym_op
+{
+ signed int *h;
+ signed int *k;
+ signed int *l; /* Contributions to h, k and l from h, k, i and l */
+ enum lattice_type latt;
+ char op[6];
};
@@ -59,416 +68,341 @@ struct sym_op {
* @include: "symmetry.h"
* @Image:
*
- * Wibble
+ * The SymOpList is an opaque data structure containing a list of point symmetry
+ * operations. It could represent an point group or a list of indexing
+ * ambiguities (twin laws), or similar.
*/
-struct _symoplist {
- struct sym_op *items;
+struct _symoplist
+{
+ struct sym_op *ops;
int n_ops;
int max_ops;
};
-static void alloc_items(SymOpList *items)
+static void alloc_ops(SymOpList *ops)
{
- items->items = realloc(items->items,
- items->max_ops*sizeof(struct sym_op));
+ ops->ops = realloc(ops->ops, ops->max_ops*sizeof(struct sym_op));
}
-/**
- * new_items:
- *
- * Creates a new %SymOpList.
- *
- * Returns: The new list, or NULL.
- **/
-SymOpList *new_items()
+/* Creates a new SymOpList */
+static SymOpList *new_symoplist()
{
SymOpList *new;
new = malloc(sizeof(SymOpList));
if ( new == NULL ) return NULL;
- new->max_ops = 1024;
+ new->max_ops = 16;
new->n_ops = 0;
- new->items = NULL;
- alloc_items(new);
+ new->ops = NULL;
+ alloc_ops(new);
return new;
}
-void delete_items(SymOpList *items)
+/**
+ * free_symoplist:
+ *
+ * Frees a %SymOpList and all associated resources.
+ **/
+void free_symoplist(SymOpList *ops)
{
- if ( items == NULL ) return;
- if ( items->items != NULL ) free(items->items);
- free(items);
+ int i;
+
+ if ( ops == NULL ) return;
+ for ( i=0; i<ops->n_ops; i++ ) {
+ free(ops->ops[i].h);
+ free(ops->ops[i].k);
+ free(ops->ops[i].l);
+ }
+ if ( ops->ops != NULL ) free(ops->ops);
+ free(ops);
}
-void add_item_with_op(SymOpList *items, signed int h, signed int k,
- signed int l, int op)
+/* Add a operation to a SymOpList */
+static void add_symop(SymOpList *ops,
+ signed int *h, signed int *k, signed int *l)
{
- if ( items->n_ops == items->max_ops ) {
- items->max_ops += 1024;
- alloc_items(items);
+ if ( ops->n_ops == ops->max_ops ) {
+ /* Pretty sure this never happens, but still... */
+ ops->max_ops += 16;
+ alloc_ops(ops);
}
- items->items[items->n_ops].h = h;
- items->items[items->n_ops].k = k;
- items->items[items->n_ops].l = l;
- items->items[items->n_ops].op = op;
- items->n_ops++;
+ ops->ops[ops->n_ops].h = h;
+ ops->ops[ops->n_ops].k = k;
+ ops->ops[ops->n_ops].l = l;
+ ops->n_ops++;
}
-void add_item(SymOpList *items, signed int h, signed int k, signed int l)
+int num_ops(const SymOpList *ops)
{
- add_item_with_op(items, h, k, l, 0);
+ return ops->n_ops;
}
-int find_item(SymOpList *items,
- signed int h, signed int k, signed int l)
+static signed int *v(signed int h, signed int k, signed int i, signed int l)
{
- int i;
+ signed int *vec = malloc(4*sizeof(signed int));
+ vec[0] = h; vec[1] = k; vec[2] = i; vec[3] = l;
+ return vec;
+}
- for ( i=0; i<items->n_ops; i++ ) {
- if ( items->items[i].h != h ) continue;
- if ( items->items[i].k != k ) continue;
- if ( items->items[i].l != l ) continue;
- return 1;
- }
- return 0;
+
+/********************************* Triclinic **********************************/
+
+static SymOpList *make_1bar()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,1));
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1));
+ return new;
}
-static int find_op(SymOpList *items, int op)
+static SymOpList *make_1()
{
- int i;
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,1));
+ return new;
+}
- for ( i=0; i<items->n_ops; i++ ) {
- if ( items->items[i].op == op ) return 1;
- }
- return 0;
+
+/********************************* Monoclinic *********************************/
+
+static SymOpList *make_2m()
+{
+ return NULL;
}
-struct sym_op *get_item(SymOpList *items, int i)
+static SymOpList *make_2()
{
- if ( i >= items->n_ops ) return NULL;
- return &items->items[i];
+ return NULL;
}
-int num_items(const SymOpList *items)
+static SymOpList *make_m()
{
- return items->n_ops;
+ return NULL;
}
-void union_op_items(SymOpList *items, SymOpList *newi)
+/******************************** Orthorhombic ********************************/
+
+static SymOpList *make_mmm()
{
- int n, i;
+ return NULL;
+}
- n = num_items(newi);
- for ( i=0; i<n; i++ ) {
- struct sym_op *r = get_item(newi, i);
- if ( find_op(items, r->op) ) continue;
+static SymOpList *make_222()
+{
+ return NULL;
+}
- add_item_with_op(items, r->h, r->k, r->l, r->op);
- }
+static SymOpList *make_mm2()
+{
+ return NULL;
}
-void union_ops(SymOpList *items, SymOpList *newi)
+/********************************* Tetragonal *********************************/
+
+static SymOpList *make_4m()
{
- int n, i;
+ return NULL;
+}
- n = num_items(newi);
- for ( i=0; i<n; i++ ) {
- struct sym_op *r = get_item(newi, i);
- if ( find_item(items, r->h, r->k, r->l) ) continue;
+static SymOpList *make_4()
+{
+ return NULL;
+}
- add_item_with_op(items, r->h, r->k, r->l, r->op);
- }
+static SymOpList *make_4bar()
+{
+ return NULL;
}
-SymOpList *intersection_ops(SymOpList *i1, SymOpList *i2)
+static SymOpList *make_4mmm()
{
- int n, i;
- SymOpList *res = new_items();
+ return NULL;
+}
- n = num_items(i1);
- for ( i=0; i<n; i++ ) {
- struct sym_op *r = get_item(i1, i);
- if ( find_item(i2, r->h, r->k, r->l) ) {
- add_item_with_op(res, r->h, r->k, r->l, r->op);
- }
+static SymOpList *make_422()
+{
+ return NULL;
+}
- }
- return res;
+static SymOpList *make_4bar2m()
+{
+ return NULL;
}
-/* Check if a reflection is in the asymmetric unit cell */
-static int check_cond(signed int h, signed int k, signed int l, const char *sym)
+static SymOpList *make_4mm()
{
- /* Triclinic */
- if ( strcmp(sym, "1") == 0 )
- return ( 1 );
- if ( strcmp(sym, "-1") == 0 )
- return ( (l>0)
- || ( (l==0) && (k>0) )
- || ( (l==0) && (k==0) && (h>=0) ) );
+ return NULL;
+}
- /* Orthorhombic */
- if ( strcmp(sym, "mmm") == 0 )
- return ( (h>=0) && (k>=0) && (l>=0) );
- if ( strcmp(sym, "222") == 0 )
- return ( (h>=0) && (k>0) )
- || ( (h>0) && (k>=0) )
- || ( (h>=0) && (k==0) && (l>=0) )
- || ( (h==0) && (k>=0) && (l>=0) );
- /* Tetragonal */
- if ( strcmp(sym, "4/mmm") == 0 )
- return ( (((h>0) && (k>=0)) || ((h==0) && (k==0))) && (l>=0)
- && (h>=k) ); /* Like 6/mmm */
- if ( strcmp(sym, "422") == 0 )
- return ( (((h>0) && (k>=0)) || ((h==0) && (k==0)))
- && (h>=k) );
- if ( strcmp(sym, "4/m") == 0 )
- return ( (((h>0) && (k>=0)) || ((h==0) && (k==0))) && (l>=0) );
- if ( strcmp(sym, "4") == 0 )
- return ( ((h>0) && (k>=0)) || ((h==0) && (k==0)) );
+/******************************** Rhombohedral ********************************/
- /* Hexagonal */
- if ( strcmp(sym, "6") == 0 )
- return ( ((h>0) && (k>=0)) || ((h==0) && (k==0)) );
- if ( strcmp(sym, "6/m") == 0 )
- return ( (((h>0) && (k>=0)) || ((h==0) && (k==0))) && (l>=0) );
- if ( strcmp(sym, "6/mmm") == 0 )
- return ( (((h>0) && (k>=0)) || ((h==0) && (k==0))) && (l>=0)
- && (h>=k) );
- /* TODO: Add more groups here */
+/********************************** Hexgonal **********************************/
- return 1;
+static SymOpList *make_6m()
+{
+ return NULL;
}
-/* Macros for checking the above conditions and returning if satisfied */
-#define CHECK_COND(h, k, l, sym) \
- if ( check_cond((h), (k), (l), (sym)) ) { \
- *hp = (h); *kp = (k); *lp = (l); \
- return; \
- }
+static SymOpList *make_6()
+{
+ return NULL;
+}
-int num_general_equivs(const char *sym)
+static SymOpList *make_6bar()
{
- /* Triclinic */
- if ( strcmp(sym, "1") == 0 ) return 1;
- if ( strcmp(sym, "-1") == 0 ) return 2;
+ return NULL;
+}
- /* Orthorhombic */
- if ( strcmp(sym, "222") == 0 ) return 4;
- if ( strcmp(sym, "mmm") == 0 ) return 8;
- /* Tetragonal */
- if ( strcmp(sym, "4") == 0 ) return 4;
- if ( strcmp(sym, "4/m") == 0 ) return 8;
- if ( strcmp(sym, "422") == 0 ) return 8;
- if ( strcmp(sym, "4/mmm") == 0 ) return 16;
+static SymOpList *make_6mmm()
+{
+ return NULL;
+}
- /* 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 */
+static SymOpList *make_622()
+{
+ return NULL;
+}
+
- return 1;
+static SymOpList *make_6bar2m()
+{
+ return NULL;
}
-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)
+static SymOpList *make_6mm()
{
- signed int i = -h-k;
+ return NULL;
+}
- /* The returned indices when idx=0 *must* be the same as the input.
- * After that, the order does not matter. */
- if ( strcmp(sym, "1") == 0 ) {
- *he = h; *ke = k; *le = l; return;
- }
+/************************************ Cubic ***********************************/
- if ( strcmp(sym, "-1") == 0 ) {
- switch ( idx ) {
- case 0 : *he = h; *ke = k; *le = l; return;
- case 1 : *he = -h; *ke = -k; *le = -l; return;
- }
- }
- if ( strcmp(sym, "222") == 0 ) {
- switch ( idx ) {
- case 0 : *he = h; *ke = k; *le = l; return;
- case 1 : *he = -h; *ke = -k; *le = l; return;
- case 2 : *he = -h; *ke = k; *le = -l; return;
- case 3 : *he = h; *ke = -k; *le = -l; return;
- }
- }
- if ( strcmp(sym, "mmm") == 0 ) {
- switch ( idx ) {
- case 0 : *he = h; *ke = k; *le = l; return;
- case 1 : *he = -h; *ke = -k; *le = l; return;
- case 2 : *he = -h; *ke = k; *le = -l; return;
- case 3 : *he = h; *ke = -k; *le = -l; return;
- case 4 : *he = -h; *ke = -k; *le = -l; return;
- case 5 : *he = h; *ke = k; *le = -l; return;
- case 6 : *he = h; *ke = -k; *le = l; return;
- case 7 : *he = -h; *ke = k; *le = l; return;
- }
- }
+SymOpList *get_pointgroup(const char *sym)
+{
+ /* Triclinic */
+ if ( strcmp(sym, "-1") == 0 ) return make_1bar();
+ if ( strcmp(sym, "1") == 0 ) return make_1();
- if ( strcmp(sym, "4") == 0 ) {
- switch ( idx ) {
- case 0 : *he = h; *ke = k; *le = l; return;
- case 1 : *he = -h; *ke = -k; *le = l; return;
- case 2 : *he = -k; *ke = h; *le = l; return;
- case 3 : *he = k; *ke = -h; *le = l; return;
- }
- }
+ /* Monoclinic */
+ if ( strcmp(sym, "2/m") == 0 ) return make_2m();
+ if ( strcmp(sym, "2") == 0 ) return make_2();
+ if ( strcmp(sym, "m") == 0 ) return make_m();
- if ( strcmp(sym, "4/m") == 0 ) {
- switch ( idx ) {
- case 0 : *he = h; *ke = k; *le = l; return;
- case 1 : *he = -h; *ke = -k; *le = l; return;
- case 2 : *he = -k; *ke = h; *le = l; return;
- case 3 : *he = k; *ke = -h; *le = l; return;
- case 4 : *he = -h; *ke = -k; *le = -l; return;
- case 5 : *he = h; *ke = k; *le = -l; return;
- case 6 : *he = k; *ke = -h; *le = -l; return;
- case 7 : *he = -k; *ke = h; *le = -l; return;
- }
- }
+ /* Orthorhombic */
+ if ( strcmp(sym, "mmm") == 0 ) return make_mmm();
+ if ( strcmp(sym, "222") == 0 ) return make_222();
+ if ( strcmp(sym, "mm2") == 0 ) return make_mm2();
- if ( strcmp(sym, "422") == 0 ) {
- switch ( idx ) {
- case 0 : *he = h; *ke = k; *le = l; return;
- case 1 : *he = -h; *ke = -k; *le = l; return;
- case 2 : *he = -k; *ke = h; *le = l; return;
- case 3 : *he = k; *ke = -h; *le = l; return;
- case 4 : *he = -h; *ke = k; *le = -l; return;
- case 5 : *he = h; *ke = -k; *le = -l; return;
- case 6 : *he = k; *ke = h; *le = -l; return;
- case 7 : *he = -k; *ke = -h; *le = -l; return;
- }
- }
+ /* Tetragonal */
+ if ( strcmp(sym, "4/m") == 0 ) return make_4m();
+ if ( strcmp(sym, "4") == 0 ) return make_4();
+ if ( strcmp(sym, "4bar") == 0 ) return make_4bar();
+ if ( strcmp(sym, "4/mmm") == 0 ) return make_4mmm();
+ if ( strcmp(sym, "422") == 0 ) return make_422();
+ if ( strcmp(sym, "4bar2m") == 0 ) return make_4bar2m();
+ if ( strcmp(sym, "4mm") == 0 ) return make_4mm();
- if ( strcmp(sym, "4/mmm") == 0 ) {
- switch ( idx ) {
- case 0 : *he = h; *ke = k; *le = l; return;
- case 1 : *he = -h; *ke = -k; *le = l; return;
- case 2 : *he = -k; *ke = h; *le = l; return;
- case 3 : *he = k; *ke = -h; *le = l; return;
- case 4 : *he = -h; *ke = k; *le = -l; return;
- case 5 : *he = h; *ke = -k; *le = -l; return;
- case 6 : *he = k; *ke = h; *le = -l; return;
- case 7 : *he = -k; *ke = -h; *le = -l; return;
- case 8 : *he = -h; *ke = -k; *le = -l; return;
- case 9 : *he = h; *ke = k; *le = -l; return;
- case 10 : *he = k; *ke = -h; *le = -l; return;
- case 11 : *he = -k; *ke = h; *le = -l; return;
- case 12 : *he = h; *ke = -k; *le = l; return;
- case 13 : *he = -h; *ke = k; *le = l; return;
- case 14 : *he = -k; *ke = -h; *le = l; return;
- case 15 : *he = k; *ke = h; *le = l; return;
- }
- }
+ /* Hexagonal */
+ if ( strcmp(sym, "6/m") == 0 ) return make_6m();
+ if ( strcmp(sym, "6") == 0 ) return make_6();
+ if ( strcmp(sym, "6bar") == 0 ) return make_6bar();
+ if ( strcmp(sym, "6/mmm") == 0 ) return make_6mmm();
+ if ( strcmp(sym, "622") == 0 ) return make_622();
+ if ( strcmp(sym, "6bar2m") == 0 ) return make_6bar2m();
+ if ( strcmp(sym, "6mm") == 0 ) return make_6mm();
+}
- if ( strcmp(sym, "6") == 0 ) {
- switch ( idx ) {
- case 0 : *he = h; *ke = k; *le = l; return;
- case 1 : *he = i; *ke = h; *le = l; return;
- case 2 : *he = k; *ke = i; *le = l; return;
- case 3 : *he = -h; *ke = -k; *le = l; return;
- case 4 : *he = -i; *ke = -h; *le = l; return;
- case 5 : *he = -k; *ke = -i; *le = l; return;
- }
- }
- if ( strcmp(sym, "6/m") == 0 ) {
- switch ( idx ) {
- case 0 : *he = h; *ke = k; *le = l; return;
- case 1 : *he = i; *ke = h; *le = l; return;
- case 2 : *he = k; *ke = i; *le = l; return;
- case 3 : *he = -h; *ke = -k; *le = l; return;
- case 4 : *he = -i; *ke = -h; *le = l; return;
- case 5 : *he = -k; *ke = -i; *le = l; return;
- case 6 : *he = h; *ke = k; *le = -l; return;
- case 7 : *he = i; *ke = h; *le = -l; return;
- case 8 : *he = k; *ke = i; *le = -l; return;
- case 9 : *he = -h; *ke = -k; *le = -l; return;
- case 10 : *he = -i; *ke = -h; *le = -l; return;
- case 11 : *he = -k; *ke = -i; *le = -l; return;
- }
- }
+/**
+ * num_ops:
+ * @ops: A %SymOpList
+ *
+ * Returns: the number of operations in @ops.
+ **/
+int num_ops(SymOpList *ops)
+{
+ return ops->n_ops;
+}
- if ( strcmp(sym, "6/mmm") == 0 ) {
- switch ( idx ) {
- case 0 : *he = h; *ke = k; *le = l; return;
- case 1 : *he = i; *ke = h; *le = l; return;
- case 2 : *he = k; *ke = i; *le = l; return;
- case 3 : *he = -h; *ke = -k; *le = l; return;
- case 4 : *he = -i; *ke = -h; *le = l; return;
- case 5 : *he = -k; *ke = -i; *le = l; return;
- case 6 : *he = k; *ke = h; *le = -l; return;
- case 7 : *he = h; *ke = i; *le = -l; return;
- case 8 : *he = i; *ke = k; *le = -l; return;
- case 9 : *he = -k; *ke = -h; *le = -l; return;
- case 10 : *he = -h; *ke = -i; *le = -l; return;
- case 11 : *he = -i; *ke = -k; *le = -l; return;
- case 12 : *he = -h; *ke = -k; *le = -l; return;
- case 13 : *he = -i; *ke = -h; *le = -l; return;
- case 14 : *he = -k; *ke = -i; *le = -l; return;
- case 15 : *he = h; *ke = k; *le = -l; return;
- case 16 : *he = i; *ke = h; *le = -l; return;
- case 17 : *he = k; *ke = i; *le = -l; return;
- case 18 : *he = -k; *ke = -h; *le = l; return;
- case 19 : *he = -h; *ke = -i; *le = l; return;
- case 20 : *he = -i; *ke = -k; *le = l; return;
- case 21 : *he = k; *ke = h; *le = l; return;
- case 22 : *he = h; *ke = i; *le = l; return;
- case 23 : *he = i; *ke = k; *le = l; return;
- }
- }
- /* TODO: Add more groups here */
+/**
+ * get_equiv:
+ * @ops: A %SymOpList
+ * @idx: Index of the operation to use
+ * @h: index of reflection
+ * @k: index of reflection
+ * @l: index of reflection
+ * @he: location to store h index of equivalent reflection
+ * @ke: location to store k index of equivalent reflection
+ * @le: location to store l index of equivalent reflection
+ *
+ * This function applies the @idx-th symmetry operation from @ops to the
+ * reflection @h, @k, @l, and stores the result at @he, @ke and @le.
+ *
+ * If you don't mind that the same equivalent might appear twice, simply call
+ * this function the number of times returned by num_ops(), using the actual
+ * point group. If repeating the same equivalent twice (for example, if the
+ * given reflection is a special high-symmetry one), call special_position()
+ * first to get a "specialised" SymOpList and use that instead.
+ **/
+void get_equiv(SymOpList *ops, int idx,
+ signed int h, signed int k, signed int l,
+ signed int *he, signed int *ke, signed int *le)
+{
+ signed int i = -h-k;
+ struct sym_op op = ops->ops[idx];
- ERROR("Unrecognised symmetry '%s'\n", sym);
- abort();
+ *he = h*op.h[0] + k*op.h[1] + i*op.h[2] + l*op.h[3];
+ *ke = h*op.k[0] + k*op.h[1] + i*op.k[2] + l*op.k[3];
+ *le = h*op.l[0] + k*op.h[1] + i*op.l[2] + l*op.l[3];
}
-/* Given a reflection and a point group, this returns (by reference) the indices
- * of the "idx"th equivalent reflection, taking special positions into account.
- * It returns "idx" if successful. Otherwise, it returns the number of
- * equivalents for the particular reflection (taking special positions into
- * account). Therefore, set idx=-1 to get the number of equivalents. */
-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)
+/**
+ * special_position:
+ * @ops: A %SymOpList, usually corresponding to a point group
+ * @h: index of a reflection
+ * @k: index of a reflection
+ * @l: index of a reflection
+ *
+ * This function determines which operations in @ops map the reflection @h, @k,
+ * @l onto itself, and returns a new %SymOpList containing only the operations
+ * from @ops which do not do so.
+ *
+ * Returns: the "specialised" %SymOpList.
+ **/
+SymOpList *special_position(SymOpList *ops,
+ signed int h, signed int k, signed int l)
{
int n_general;
int i;
@@ -481,15 +415,15 @@ static int special_position(signed int hs, signed int ks, signed int ls,
return 0;
}
- equivs = new_items();
- n_general = num_general_equivs(sym);
+ equivs = new_symoplist();
+ n_general = num_general_equivs(ops);
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);
+ /* Get equivalent according to the point group */
+ get_general_equiv(ops, i, hs, ks, ls, &h, &k, &l);
/* Already got this one? */
if ( find_item(equivs, h, k, l) ) continue;
@@ -498,7 +432,7 @@ static int special_position(signed int hs, signed int ks, signed int ls,
*hp = h;
*kp = k;
*lp = l;
- delete_items(equivs);
+ delete_ops(equivs);
return n_equivs;
}
add_item(equivs, h, k, l);
@@ -506,329 +440,47 @@ static int special_position(signed int hs, signed int ks, signed int ls,
}
- delete_items(equivs);
+ delete_ops(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)
+void get_asymm(SymOpList *ops, int idx,
+ signed int h, signed int k, signed int l,
+ signed int *hp, signed int *kp, signed int *lp)
{
int nequiv = num_equivs(h, k, l, sym);
int p;
- SYM_DEBUG("------ %i %i %i\n", h, k, l);
for ( p=0; p<nequiv; p++ ) {
signed int he, ke, le;
get_equiv(h, k, l, &he, &ke, &le, sym, p);
- SYM_DEBUG("%i : %i %i %i\n", p, he, ke, le);
CHECK_COND(he, ke, le, sym);
}
-
- /* Should never reach here */
- ERROR("No match found in %s for %i %i %i\n", sym, h, k, l);
- abort();
}
-/* This is kind of like a "numerical" left coset decomposition.
- * Given a reflection index and a point group, it returns the list of twinning
- * possibilities.
+/**
+ * get_twins:
*
- * To count the number of possibilities, use num_items() on the result.
+ * Calculate twinning laws.
+ *
+ * To count the number of possibilities, use num_ops() on the result.
*/
-static SymOpList *coset_decomp(signed int hs, signed int ks, signed int ls,
- const char *holo, const char *mero)
+SymOpList *get_twins(SymOpList *source, SymOpList *target)
{
- int n_mero, n_holo;
+ int n_src, n_tgt;
int i;
signed int h, k, l;
- SymOpList *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_general_equivs(holo);
- n_mero = num_general_equivs(mero);
-
- for ( i=0; i<n_holo; i++ ) {
+ SymOpList *twins = new_symoplist();
- signed int h_holo, k_holo, l_holo;
- signed int hs_holo, ks_holo, ls_holo;
+ n_src = num_ops(source);
+ n_tgt = num_ops(target);
- /* Get equivalent according to the holohedral group */
- get_general_equiv(h, k, l, &hs_holo, &ks_holo, &ls_holo,
- holo, i);
+ for ( i=0; i<n_src; 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);
-
- /* 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);
}
return twins;
}
-
-
-/* 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. */
-SymOpList *get_twins(const char *holo, const char *mero)
-{
- int i;
- SymOpList *ops = new_items();
- int expected, actual;
- SymOpList *items;
-
- /* 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 sym_op *item;
- SymOpList *new_ops;
-
- item = get_item(items, i);
-
- h = item->h;
- k = item->k;
- l = item->l;
-
- new_ops = coset_decomp(h, k, l, holo, mero);
- union_op_items(ops, new_ops);
- delete_items(new_ops);
-
- }
-
- /* Idiot check */
- actual = num_items(ops);
- expected = num_general_equivs(holo) / num_general_equivs(mero);
- if ( actual != expected ) {
- ERROR("Found %i twin possibilities, but expected %i.\n",
- actual, expected);
- ERROR("I couldn't find the number of twin laws that I expected."
- " This is an internal error, and shouldn't happen. "
- "Sorry.\n");
- abort();
- }
-
- return ops;
-}
-
-
-int find_unique_equiv(SymOpList *items, signed int h, signed int k,
- signed int l, const char *mero, signed int *hu,
- signed int *ku, signed int *lu)
-{
- int i;
- int found = 0;
-
- for ( i=0; i<num_equivs(h, k, l, mero); i++ ) {
-
- signed int he, ke, le;
- int f;
- get_equiv(h, k, l, &he, &ke, &le, mero, i);
- 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.
- * 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;
- return 1;
- }
-
- }
-
- return 0;
-}
-
-
-/* Returns true if the point group is 23, m-3, 432, -43m or m-3m
- * (i.e. T, Th, O, Td or Oh) */
-int is_polyhedral(const char *sym)
-{
- /* Triclinic */
- if ( strcmp(sym, "1") == 0 ) return 0;
- if ( strcmp(sym, "-1") == 0 ) return 0;
-
- /* Orthorhombic */
- if ( strcmp(sym, "222") == 0 ) return 0;
- if ( strcmp(sym, "mmm") == 0 ) return 0;
-
- /* Tetragonal */
- if ( strcmp(sym, "4") == 0 ) return 0;
- if ( strcmp(sym, "4/m") == 0 ) return 0;
- if ( strcmp(sym, "422") == 0 ) return 0;
- if ( strcmp(sym, "4/mmm") == 0 ) return 0;
-
- /* Hexagonal */
- if ( strcmp(sym, "6") == 0 ) return 0;
- if ( strcmp(sym, "6/m") == 0 ) return 0;
- if ( strcmp(sym, "6/mmm") == 0 ) return 0;
-
- /* TODO: Add more groups here */
-
- ERROR("Don't know if '%s' is polyhedral or not.\n", sym);
- abort();
-}
-
-
-/* Returns the order of the characteristic axis of proper or improper rotation
- * for the point group. This should be the rotation about the c axis. */
-int rotational_order(const char *sym)
-{
- /* Triclinic */
- if ( strcmp(sym, "1") == 0 ) return 1;
- if ( strcmp(sym, "-1") == 0 ) return 2;
-
- /* Orthorhombic */
- if ( strcmp(sym, "222") == 0 ) return 2;
- if ( strcmp(sym, "mmm") == 0 ) return 2;
-
- /* Tetragonal */
- if ( strcmp(sym, "4") == 0 ) return 4;
- if ( strcmp(sym, "4/m") == 0 ) return 4;
- if ( strcmp(sym, "422") == 0 ) return 4;
- if ( strcmp(sym, "4/mmm") == 0 ) return 4;
-
- /* Hexagonal */
- if ( strcmp(sym, "6") == 0 ) return 6;
- if ( strcmp(sym, "6/m") == 0 ) return 6;
- if ( strcmp(sym, "6/mmm") == 0 ) return 6;
-
- /* TODO: Add more groups here */
-
- ERROR("Couldn't find rotational order for '%s'.\n", sym);
- abort();
-}
-
-
-int has_perpendicular_mirror(const char *sym)
-{
- /* Triclinic */
- if ( strcmp(sym, "1") == 0 ) return 0;
- if ( strcmp(sym, "-1") == 0 ) return 0;
-
- /* Orthorhombic */
- if ( strcmp(sym, "222") == 0 ) return 0;
- if ( strcmp(sym, "mmm") == 0 ) return 1;
-
- /* Tetragonal */
- if ( strcmp(sym, "4") == 0 ) return 0;
- if ( strcmp(sym, "4/m") == 0 ) return 1;
- if ( strcmp(sym, "422") == 0 ) return 0;
- if ( strcmp(sym, "4/mmm") == 0 ) return 1;
-
- /* Hexagonal */
- if ( strcmp(sym, "6") == 0 ) return 0;
- if ( strcmp(sym, "6/m") == 0 ) return 1;
- if ( strcmp(sym, "6/mmm") == 0 ) return 1;
-
- /* TODO: Add more groups here */
-
- ERROR("Couldn't find mirror definition for '%s'.\n", sym);
- abort();
-}
-
-
-int has_bisecting_mirror_or_diad(const char *sym)
-{
- /* Triclinic */
- if ( strcmp(sym, "1") == 0 ) return 0;
- if ( strcmp(sym, "-1") == 0 ) return 0;
-
- /* Orthorhombic */
- if ( strcmp(sym, "222") == 0 ) return 1;
- if ( strcmp(sym, "mmm") == 0 ) return 1;
-
- /* Tetragonal */
- if ( strcmp(sym, "4") == 0 ) return 0;
- if ( strcmp(sym, "4/m") == 0 ) return 0;
- if ( strcmp(sym, "422") == 0 ) return 0;
- if ( strcmp(sym, "4/mmm") == 0 ) return 1;
-
- /* Hexagonal */
- if ( strcmp(sym, "6") == 0 ) return 0;
- if ( strcmp(sym, "6/m") == 0 ) return 1;
- if ( strcmp(sym, "6/mmm") == 0 ) return 1;
-
- /* TODO: Add more groups here */
-
- ERROR("Couldn't find mirror definition for '%s'.\n", sym);
- abort();
-}
-
-
-int check_symmetry(SymOpList *items, const char *sym)
-{
- int i;
- unsigned char *flags;
-
- flags = new_list_flag();
- for ( i=0; i<num_items(items); i++ ) {
- struct sym_op *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 sym_op *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 ) {
- free(flags);
- return 1; /* Symmetry is wrong! */
- }
-
- }
-
- free(flags);
-
- return 0;
-}
diff --git a/src/symmetry.h b/src/symmetry.h
index fbb80d26..23c00160 100644
--- a/src/symmetry.h
+++ b/src/symmetry.h
@@ -24,6 +24,12 @@
**/
typedef struct _symoplist SymOpList;
+extern void free_symoplist(SymOpList *ops);
+
+extern SymOpList *get_pointgroup(const char *sym);
+
+extern SymOpList *special_position(SymOpList *ops,
+ signed int h, signed int k, signed int l);
extern void get_asymm(signed int h, signed int k, signed int l,
signed int *hp, signed int *kp, signed int *lp,
@@ -42,12 +48,6 @@ 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 SymOpList *get_twins(const char *holo, const char *mero);
-
-/* Properties of point groups */
-extern int is_polyhedral(const char *sym);
-extern int rotational_order(const char *sym);
-extern int has_perpendicular_mirror(const char *sym);
-extern int has_bisecting_mirror_or_diad(const char *sym);
+extern SymOpList *get_twins(SymOpList *source, SymOpList *target);
#endif /* SYMMETRY_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-04 15:39:54 +0000