aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--doc/symmetry59
-rw-r--r--src/symmetry.c229
-rw-r--r--src/symmetry.h10
-rw-r--r--src/utils.c21
-rw-r--r--src/utils.h3
5 files changed, 176 insertions, 146 deletions
diff --git a/doc/symmetry b/doc/symmetry
index 61bc98e3..d201e000 100644
--- a/doc/symmetry
+++ b/doc/symmetry
@@ -15,14 +15,15 @@ model using something like:
$ echo HYDROGENS APPEND | hgen xyzin model.pdb xyzout model-with-H.pdb ]
Symmetry definitions are included in src/symmetry.c. Point group definitions
-are required for merging and the display of merged results (since merging does
-not care about systematic absences (i.e. the space group) - as far as merging is
-concerned, systematic absences are just measurements which happen to have values
-of zero. Each space group belongs to exactly one point group, which you can
-look up in the International Tables for X-Ray Crystallography.
-
-Space groups would be needed to make get_hkl handle symmetry properly, but that
-hasn't been done yet, so symmetry.c just handles point groups for now. The
+are required for merging and the display of merged results, but space groups are
+not taken into account since merging does not care about systematic absences -
+as far as process_hkl is concerned, systematic absences are just measurements
+which happen to have values of zero. Each space group belongs to exactly one
+point group, which you can look up in the International Tables for X-Ray
+Crystallography.
+
+Space groups would only be needed to make get_hkl handle symmetry properly, but
+that hasn't been done yet, so symmetry.c just handles point groups for now. The
method used in symmetry.c is general to both point groups and space groups, even
though the code currently is not.
@@ -35,30 +36,18 @@ whole of reciprocal space must map onto exactly one reflection in the asymmetric
unit cell so defined. The asymmetric cell is usually defined with positive h, k
and l, but it doesn't really matter. Working out the required condition means
visualising the cell and taking care to properly handle situations such as the
-(000) reflection.
-
-Now, expand the num_equivs() function. Given a point group and a set of Miller
-indices, this function must return the number of equivalent reflections,
-including the given reflection, taking into account the multiplicity of the
-reflection. For example, high-symmetry reflections (usually ones with zeroes
-in their indices) have fewer equivalents. Label the different classes of
-reflection according to their Wyckoff letters. This information can be found in
-the International Tables.
-
-Finally, add the new point group to the get_equiv() function. This function
-takes a set of Miller indices, a point group and an index, "n", and returns (by
-reference) the indices of the "n"th equivalent reflection. Again, you have to
-distinguish between the different Wyckoff positions. Lists of equivalent
-reflection, grouped according to Wyckoff symbol, can also be found in the
-International Tables. get_equiv() should return the original indices as the
-first set of equivalent indices, i.e. when idx=0.
-
-It's not normally necessary to write an individual table of equivalents for each
-Wyckoff position, because many positions are just "truncated" sub-classes of
-other positions. For example, in 6/mmm, the equivalent reflections for "b" and
-"c" can be generated simply by working through the list for "g" (the general
-position) with indices set to zero or equal as appropriate, and by stopping
-after the first six equivalents. Therefore, the num_equivs() function combined
-with a single table for positions b-g in get_equivs() is sufficient to deal with
-all cases. Check carefully whether this really works for your chosen point
-group first.
+(000) reflection. Get this right, otherwise you'll go crazy when your symmetry
+breaks in weird ways.
+
+Now, expand the num_general_equivs() function. Given a point group, this
+function must return the number of equivalent reflections for a general
+reflection, including the given reflection. High-symmetry reflections (usually
+ones with zeroes in their indices) have fewer equivalents, but the num_equivs()
+function will work this out for you.
+
+Finally, add the new point group to the get_general_equiv() function. This
+function takes a set of Miller indices, a point group and an index "n", and
+returns (by reference) the indices of the "n"th equivalent reflection. You just
+have to worry about the general position, because get_equiv() will work out the
+special positions for you. get_general_equiv() should return the original
+indices when idx=0.
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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-02 13:55:22 +0000