aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-07-20 17:53:27 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:33 +0100
commitb143429764665a75dd3baf8c5115bf8553d18d71 (patch)
treeaee21b359e82c9c107e5f82de43c2d703c45538b
parent012073a3be1bb523588b83d8be0589a5d00676aa (diff)
Symmetry stuff
-rw-r--r--Makefile.am9
-rw-r--r--Makefile.in28
-rw-r--r--src/diffraction.c8
-rw-r--r--src/get_hkl.c29
-rw-r--r--src/povray.c22
-rw-r--r--src/powder_plot.c10
-rw-r--r--src/process_hkl.c3
-rw-r--r--src/reflist-utils.c23
-rw-r--r--src/render_hkl.c14
-rw-r--r--src/symmetry.c423
-rw-r--r--src/symmetry.h21
11 files changed, 382 insertions, 208 deletions
diff --git a/Makefile.am b/Makefile.am
index a2fdb3c0..f9615e37 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -9,11 +9,12 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \
contrib/alter_stream
noinst_PROGRAMS = tests/list_check tests/integration_check \
- tests/pr_gradient_check
+ tests/pr_gradient_check tests/symmetry_check
TESTS = tests/list_check tests/first_merge_check tests/second_merge_check \
tests/third_merge_check tests/fourth_merge_check \
- tests/integration_check tests/pr_gradient_check
+ tests/integration_check tests/pr_gradient_check \
+ tests/symmetry_check
EXTRA_DIST += tests/first_merge_check tests/second_merge_check \
tests/third_merge_check tests/fourth_merge_check \
@@ -146,6 +147,10 @@ tests_integration_check_SOURCES = tests/integration_check.c src/peaks.c \
src/thread-pool.c src/utils.c src/image.c \
src/hdf5-file.c
+tests_symmetry_check_SOURCES = tests/symmetry_check.c src/symmetry.c \
+ src/utils.c src/thread-pool.c
+
+
tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c src/detector.c \
src/cell.c src/geometry.c src/reflist.c \
src/thread-pool.c src/utils.c src/peaks.c \
diff --git a/Makefile.in b/Makefile.in
index df5dd25d..b94918e7 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -44,11 +44,13 @@ bin_PROGRAMS = src/pattern_sim$(EXEEXT) src/process_hkl$(EXEEXT) \
contrib/alter_stream$(EXEEXT) $(am__EXEEXT_1) $(am__EXEEXT_2)
noinst_PROGRAMS = tests/list_check$(EXEEXT) \
tests/integration_check$(EXEEXT) \
- tests/pr_gradient_check$(EXEEXT) $(am__EXEEXT_3)
+ tests/pr_gradient_check$(EXEEXT) tests/symmetry_check$(EXEEXT) \
+ $(am__EXEEXT_3)
TESTS = tests/list_check$(EXEEXT) tests/first_merge_check \
tests/second_merge_check tests/third_merge_check \
tests/fourth_merge_check tests/integration_check$(EXEEXT) \
- tests/pr_gradient_check$(EXEEXT) $(am__EXEEXT_3)
+ tests/pr_gradient_check$(EXEEXT) tests/symmetry_check$(EXEEXT) \
+ $(am__EXEEXT_3)
@BUILD_HDFSEE_TRUE@am__append_1 = src/hdfsee
@BUILD_CUBEIT_TRUE@am__append_2 = src/cubeit
@HAVE_OPENCL_TRUE@am__append_3 = src/diffraction-gpu.c src/cl-utils.c
@@ -306,6 +308,12 @@ tests_pr_gradient_check_OBJECTS = \
$(am_tests_pr_gradient_check_OBJECTS)
tests_pr_gradient_check_LDADD = $(LDADD)
tests_pr_gradient_check_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
+am_tests_symmetry_check_OBJECTS = tests/symmetry_check.$(OBJEXT) \
+ src/symmetry.$(OBJEXT) src/utils.$(OBJEXT) \
+ src/thread-pool.$(OBJEXT)
+tests_symmetry_check_OBJECTS = $(am_tests_symmetry_check_OBJECTS)
+tests_symmetry_check_LDADD = $(LDADD)
+tests_symmetry_check_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
DEFAULT_INCLUDES = -I.@am__isrc@
depcomp = $(SHELL) $(top_srcdir)/build-aux/depcomp
am__depfiles_maybe = depfiles
@@ -336,7 +344,8 @@ SOURCES = $(contrib_alter_stream_SOURCES) \
$(src_render_hkl_SOURCES) $(src_sum_stack_SOURCES) \
$(tests_gpu_sim_check_SOURCES) \
$(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \
- $(tests_pr_gradient_check_SOURCES)
+ $(tests_pr_gradient_check_SOURCES) \
+ $(tests_symmetry_check_SOURCES)
DIST_SOURCES = $(contrib_alter_stream_SOURCES) \
$(src_calibrate_detector_SOURCES) $(src_check_hkl_SOURCES) \
$(src_compare_hkl_SOURCES) $(am__src_cubeit_SOURCES_DIST) \
@@ -348,7 +357,8 @@ DIST_SOURCES = $(contrib_alter_stream_SOURCES) \
$(src_sum_stack_SOURCES) \
$(am__tests_gpu_sim_check_SOURCES_DIST) \
$(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \
- $(tests_pr_gradient_check_SOURCES)
+ $(tests_pr_gradient_check_SOURCES) \
+ $(tests_symmetry_check_SOURCES)
RECURSIVE_TARGETS = all-recursive check-recursive dvi-recursive \
html-recursive info-recursive install-data-recursive \
install-dvi-recursive install-exec-recursive \
@@ -751,6 +761,9 @@ tests_integration_check_SOURCES = tests/integration_check.c src/peaks.c \
src/thread-pool.c src/utils.c src/image.c \
src/hdf5-file.c
+tests_symmetry_check_SOURCES = tests/symmetry_check.c src/symmetry.c \
+ src/utils.c src/thread-pool.c
+
tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c src/detector.c \
src/cell.c src/geometry.c src/reflist.c \
src/thread-pool.c src/utils.c src/peaks.c \
@@ -1040,6 +1053,11 @@ tests/pr_gradient_check.$(OBJEXT): tests/$(am__dirstamp) \
tests/pr_gradient_check$(EXEEXT): $(tests_pr_gradient_check_OBJECTS) $(tests_pr_gradient_check_DEPENDENCIES) tests/$(am__dirstamp)
@rm -f tests/pr_gradient_check$(EXEEXT)
$(AM_V_CCLD)$(LINK) $(tests_pr_gradient_check_OBJECTS) $(tests_pr_gradient_check_LDADD) $(LIBS)
+tests/symmetry_check.$(OBJEXT): tests/$(am__dirstamp) \
+ tests/$(DEPDIR)/$(am__dirstamp)
+tests/symmetry_check$(EXEEXT): $(tests_symmetry_check_OBJECTS) $(tests_symmetry_check_DEPENDENCIES) tests/$(am__dirstamp)
+ @rm -f tests/symmetry_check$(EXEEXT)
+ $(AM_V_CCLD)$(LINK) $(tests_symmetry_check_OBJECTS) $(tests_symmetry_check_LDADD) $(LIBS)
mostlyclean-compile:
-rm -f *.$(OBJEXT)
@@ -1089,6 +1107,7 @@ mostlyclean-compile:
-rm -f tests/integration_check.$(OBJEXT)
-rm -f tests/list_check.$(OBJEXT)
-rm -f tests/pr_gradient_check.$(OBJEXT)
+ -rm -f tests/symmetry_check.$(OBJEXT)
distclean-compile:
-rm -f *.tab.c
@@ -1139,6 +1158,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/integration_check.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/list_check.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/pr_gradient_check.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/symmetry_check.Po@am__quote@
.c.o:
@am__fastdepCC_TRUE@ $(AM_V_CC)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.o$$||'`;\
diff --git a/src/diffraction.c b/src/diffraction.c
index aa04fd52..4133cd0b 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -100,14 +100,14 @@ static double sym_lookup_intensity(const double *intensities,
int i;
double ret = 0.0;
- for ( i=0; i<num_equivs(sym); i++ ) {
+ for ( i=0; i<num_equivs(sym, NULL); i++ ) {
signed int he;
signed int ke;
signed int le;
double f, val;
- get_equiv(sym, i, h, k, l, &he, &ke, &le);
+ get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
f = (double)lookup_flag(flags, he, ke, le);
val = lookup_intensity(intensities, he, ke, le);
@@ -127,14 +127,14 @@ static double sym_lookup_phase(const double *phases,
int i;
double ret = 0.0;
- for ( i=0; i<num_equivs(sym); i++ ) {
+ for ( i=0; i<num_equivs(sym, NULL); i++ ) {
signed int he;
signed int ke;
signed int le;
double f, val;
- get_equiv(sym, i, h, k, l, &he, &ke, &le);
+ get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
f = (double)lookup_flag(flags, he, ke, le);
val = lookup_phase(phases, he, ke, le);
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 72e1c2a3..5edfa088 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -146,15 +146,17 @@ static RefList *twin_reflections(RefList *in,
Reflection *refl;
RefListIterator *iter;
RefList *out;
+ SymOpMask *m;
/* FIXME: Check properly by coset decomposition */
- if ( num_equivs(holo) < num_equivs(mero) ) {
+ if ( num_equivs(holo, NULL) < num_equivs(mero, NULL) ) {
ERROR("%s is not a subgroup of %s!\n", symmetry_name(mero),
symmetry_name(holo));
return NULL;
}
out = reflist_new();
+ m = new_symopmask(holo);
for ( refl = first_refl(in, &iter);
refl != NULL;
@@ -178,13 +180,15 @@ static RefList *twin_reflections(RefList *in,
total = 0.0;
sigma = 0.0;
skip = 0;
- n = num_equivs(holo);
+ special_position(holo, m, h, k, l);
+ n = num_equivs(holo, m);
+
for ( j=0; j<n; j++ ) {
signed int he, ke, le;
signed int hu, ku, lu;
- get_equiv(holo, j, h, k, l, &he, &ke, &le);
+ get_equiv(holo, m, j, h, k, l, &he, &ke, &le);
/* Do we have this reflection?
* We might not have the particular (merohedral)
@@ -229,15 +233,17 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
Reflection *refl;
RefListIterator *iter;
RefList *out;
+ SymOpMask *m;
/* FIXME: Check properly */
- if ( num_equivs(target) > num_equivs(initial) ) {
+ if ( num_equivs(target, NULL) > num_equivs(initial, NULL) ) {
ERROR("%s is not a subgroup of %s!\n", symmetry_name(initial),
symmetry_name(target));
return NULL;
}
out = reflist_new();
+ m = new_symopmask(initial);
for ( refl = first_refl(in, &iter);
refl != NULL;
@@ -250,7 +256,8 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
get_indices(refl, &h, &k, &l);
intensity = get_intensity(refl);
- n = num_equivs(initial);
+ special_position(initial, m, h, k, l);
+ n = num_equivs(initial, m);
/* For each equivalent in the higher symmetry group */
for ( j=0; j<n; j++ ) {
@@ -259,7 +266,7 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
Reflection *new;
/* Get the equivalent */
- get_equiv(initial, j, h, k, l, &he, &ke, &le);
+ get_equiv(initial, m, j, h, k, l, &he, &ke, &le);
/* Put it into the asymmetric unit for the target */
get_asymm(target, he, ke, le, &he, &ke, &le);
@@ -272,6 +279,8 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
}
+ free_symopmask(m);
+
return out;
}
@@ -460,6 +469,9 @@ int main(int argc, char *argv[])
Reflection *refl;
RefListIterator *iter;
+ SymOpMask *m;
+
+ m = new_symopmask(mero);
for ( refl = first_refl(input, &iter);
refl != NULL;
@@ -471,10 +483,13 @@ int main(int argc, char *argv[])
get_indices(refl, &h, &k, &l);
inty = get_intensity(refl);
- inty *= (double)num_equivs(mero);
+ special_position(mero, m, h, k, l);
+ inty *= (double)num_equivs(mero, m);
set_int(refl, inty);
}
+
+ free_symopmask(m);
}
if ( template ) {
diff --git a/src/povray.c b/src/povray.c
index 9fb267f4..d20c18e5 100644
--- a/src/povray.c
+++ b/src/povray.c
@@ -42,6 +42,7 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
int i;
Reflection *refl;
RefListIterator *iter;
+ SymOpMask *m;
if ( (nproc > MAX_PROC) || (nproc < 1) ) {
ERROR("Number of processes must be a number between 1 and %i\n",
@@ -167,6 +168,8 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
fprintf(fh, "}\n");
+ m = new_symopmask(sym);
+
max = 0.0;
for ( refl = first_refl(list, &iter);
refl != NULL;
@@ -174,10 +177,9 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
float val;
signed int h, k, l;
- SymOpList *sp;
get_indices(refl, &h, &k, &l);
- sp = special_position(sym, h, k, l);
+ special_position(sym, m, h, k, l);
switch ( wght ) {
case WGHT_I :
@@ -189,7 +191,7 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
break;
case WGHT_COUNTS :
val = get_redundancy(refl);
- val /= (double)num_equivs(sp);
+ val /= (double)num_equivs(sym, m);
break;
case WGHT_RAWCOUNTS :
val = get_redundancy(refl);
@@ -201,8 +203,6 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
if ( val > max ) max = val;
- free_symoplist(sp);
-
}
max /= boost;
@@ -211,7 +211,6 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
max = scale_top;
}
-
for ( refl = first_refl(list, &iter);
refl != NULL;
refl = next_refl(refl, iter) ) {
@@ -220,12 +219,11 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
int s;
float val, p, r, g, b, trans;
int j;
- SymOpList *sp;
int neq;
get_indices(refl, &h, &k, &l);
- sp = special_position(sym, h, k, l);
- neq = num_equivs(sp);
+ special_position(sym, m, h, k, l);
+ neq = num_equivs(sym, m);
switch ( wght ) {
case WGHT_I :
@@ -294,7 +292,7 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
signed int he, ke, le;
float x, y, z;
- get_equiv(sp, j, h, k, l, &he, &ke, &le);
+ get_equiv(sym, m, j, h, k, l, &he, &ke, &le);
x = asx*he + bsx*ke + csx*le;
y = asy*he + bsy*ke + csy*le;
@@ -311,10 +309,10 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
}
- free_symoplist(sp);
-
}
+ free_symopmask(m);
+
fprintf(fh, "\n");
fclose(fh);
diff --git a/src/powder_plot.c b/src/powder_plot.c
index 6065da07..d073f2ef 100644
--- a/src/powder_plot.c
+++ b/src/powder_plot.c
@@ -318,6 +318,9 @@ static unsigned int process_hkl(struct image *image, const SymOpList *sym,
int h, k, l, redundancy;
double q, intensity;
unsigned int nref;
+ SymOpMask *m;
+
+ m = new_symopmask(sym);
nref = num_reflections(image->reflections);
@@ -330,9 +333,8 @@ static unsigned int process_hkl(struct image *image, const SymOpList *sym,
if ( use_redundancy ) {
redundancy = get_redundancy(refl);
} else {
- SymOpList *sp = special_position(sym, h, k, l);
- redundancy = num_equivs(sp);
- free_symoplist(sp);
+ special_position(sym, m, h, k, l);
+ redundancy = num_equivs(sym, m);
}
/* Multiply by 2 to get 1/d (in m^-1) */
@@ -348,6 +350,8 @@ static unsigned int process_hkl(struct image *image, const SymOpList *sym,
}
+ free_symopmask(m);
+
return n_peaks;
}
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 06c5a54b..d34ea60a 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -521,6 +521,7 @@ int main(int argc, char *argv[])
if ( sym_str == NULL ) sym_str = strdup("1");
sym = get_pointgroup(sym_str);
+ STATUS("%s -> %p\n", sym_str, sym);
free(sym_str);
/* Open the data stream */
@@ -553,7 +554,7 @@ int main(int argc, char *argv[])
ERROR("Invalid indices for '--histogram'\n");
return 1;
}
- space_for_hist = n_total_patterns * num_equivs(sym);
+ space_for_hist = n_total_patterns * num_equivs(sym, NULL);
hist_vals = malloc(space_for_hist * sizeof(double));
free(histo);
STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l);
diff --git a/src/reflist-utils.c b/src/reflist-utils.c
index c0784dbb..5a9ec44f 100644
--- a/src/reflist-utils.c
+++ b/src/reflist-utils.c
@@ -109,8 +109,14 @@ int check_list_symmetry(RefList *list, const SymOpList *sym)
unsigned char *flags;
Reflection *refl;
RefListIterator *iter;
+ SymOpMask *mask;
flags = flags_from_list(list);
+ mask = new_symopmask(sym);
+ if ( mask == NULL ) {
+ ERROR("Couldn't create mask for list symmetry check.\n");
+ return 1;
+ }
for ( refl = first_refl(list, &iter);
refl != NULL;
@@ -119,13 +125,19 @@ int check_list_symmetry(RefList *list, const SymOpList *sym)
int j;
int found = 0;
signed int h, k, l;
+ int n;
get_indices(refl, &h, &k, &l);
- for ( j=0; j<num_equivs(sym); j++ ) {
+ special_position(sym, mask, h, k, l);
+ n = num_equivs(sym, mask);
+ STATUS("%i equivs: %3i %3i %3i\n", n, h, k, l);
+
+ for ( j=0; j<n; j++ ) {
signed int he, ke, le;
- get_equiv(sym, j, h, k, l, &he, &ke, &le);
+ get_equiv(sym, mask, j, h, k, l, &he, &ke, &le);
+ STATUS("%3i: %3i %3i %3i\n", j, he, ke, le);
if ( abs(he) > INDMAX ) continue;
if ( abs(le) > INDMAX ) continue;
@@ -137,12 +149,15 @@ int check_list_symmetry(RefList *list, const SymOpList *sym)
if ( found > 1 ) {
free(flags);
+ free_symopmask(mask);
+ STATUS("Found %i %i %i twice\n", h, k, l);
return 1; /* Symmetry is wrong! */
}
}
free(flags);
+ free_symopmask(mask);
return 0;
}
@@ -155,11 +170,11 @@ int find_equiv_in_list(RefList *list, signed int h, signed int k,
int i;
int found = 0;
- for ( i=0; i<num_equivs( sym); i++ ) {
+ for ( i=0; i<num_equivs(sym, NULL); i++ ) {
signed int he, ke, le;
Reflection *f;
- get_equiv(sym, i, h, k, l, &he, &ke, &le);
+ get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
f = find_refl(list, he, ke, le);
/* There must only be one equivalent. If there are more, it
diff --git a/src/render_hkl.c b/src/render_hkl.c
index 823b6f1a..18d6bb81 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -96,12 +96,15 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
{
Reflection *refl;
RefListIterator *iter;
+ SymOpMask *m;
if ( dctx == NULL ) {
*max_u = 0.0; *max_v = 0.0; *max_val = 0.0;
*max_res = 0.0; *max_ux = 0; *max_uy = 0;
}
+ m = new_symopmask(sym);
+
/* Iterate over all reflections */
for ( refl = first_refl(list, &iter);
refl != NULL;
@@ -111,18 +114,17 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
signed int ha, ka, la;
int xi, yi;
int i, n;
- SymOpList *sp;
get_indices(refl, &ha, &ka, &la);
- sp = special_position(sym, ha, ka, la);
- n = num_equivs(sp);
+ special_position(sym, m, ha, ka, la);
+ n = num_equivs(sym, m);
for ( i=0; i<n; i++ ) {
signed int h, k, l;
- get_equiv(sp, i, ha, ka, la, &h, &k, &l);
+ get_equiv(sym, m, i, ha, ka, la, &h, &k, &l);
/* Is the reflection in the zone? */
if ( h*zh + k*zk + l*zl != 0 ) continue;
@@ -194,9 +196,9 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
}
- free_symoplist(sp);
-
}
+
+ free_symopmask(m);
}
diff --git a/src/symmetry.c b/src/symmetry.c
index 19daf0ec..3b995967 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -37,19 +37,6 @@
*/
-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;
@@ -84,6 +71,13 @@ struct _symoplist
};
+struct _symopmask
+{
+ const SymOpList *list;
+ int *mask;
+};
+
+
static void alloc_ops(SymOpList *ops)
{
@@ -92,6 +86,35 @@ static void alloc_ops(SymOpList *ops)
}
+/**
+ * new_symopmask:
+ *
+ * Returns: a new %SymOpMask, which you can use when filtering out special
+ * reflections.
+ **/
+SymOpMask *new_symopmask(const SymOpList *list)
+{
+ SymOpMask *m;
+ int i;
+
+ m = malloc(sizeof(struct _symopmask));
+ if ( m == NULL ) return NULL;
+
+ m->list = list;
+ m->mask = malloc(sizeof(int)*list->n_ops);
+ if ( m->mask == NULL ) {
+ free(m);
+ return NULL;
+ }
+
+ for ( i=0; i<list->n_ops; i++ ) {
+ m->mask[i] = 1;
+ }
+
+ return m;
+}
+
+
/* Creates a new SymOpList */
static SymOpList *new_symoplist()
{
@@ -101,6 +124,7 @@ static SymOpList *new_symoplist()
new->max_ops = 16;
new->n_ops = 0;
new->ops = NULL;
+ new->divisors = NULL;
new->name = NULL;
new->num_equivs = 1;
alloc_ops(new);
@@ -128,88 +152,33 @@ void free_symoplist(SymOpList *ops)
free(ops);
}
-
-static int is_identity(signed int *h, signed int *k, signed int *l)
+/**
+ * free_symopmask:
+ *
+ * Frees a %SymOpMask and all associated resources.
+ **/
+void free_symopmask(SymOpMask *m)
{
- if ( (h[0]!=1) || (h[1]!=0) || (h[2]!=0) ) return 0;
- if ( (k[0]!=0) || (k[1]!=1) || (k[2]!=0) ) return 0;
- if ( (l[0]!=0) || (l[1]!=0) || (l[2]!=1) ) return 0;
- return 1;
+ if ( m == NULL ) return;
+ free(m->mask);
+ free(m);
}
-/* Calculate the order of the operation "M", which is the lowest
- * integer n such that M^n = I. */
-static int order_of_op(signed int *hin, signed int *kin, signed int *lin)
-{
- int n;
- signed int h[3];
- signed int k[3];
- signed int l[3];
-
- memcpy(h, hin, 3*sizeof(signed int));
- memcpy(k, kin, 3*sizeof(signed int));
- memcpy(l, lin, 3*sizeof(signed int));
-
- for ( n=1; n<6; n++ ) {
-
- signed int hnew[3];
- signed int knew[3];
- signed int lnew[3];
-
- /* Yay matrices */
- hnew[0] = h[0]*h[0] + h[1]*k[0] + h[2]*l[0];
- hnew[1] = h[0]*h[1] + h[1]*k[1] + h[2]*l[1];
- hnew[2] = h[0]*h[2] + h[1]*k[2] + h[2]*l[2];
-
- knew[0] = k[0]*h[0] + k[1]*k[0] + k[2]*l[0];
- knew[1] = k[0]*h[1] + k[1]*k[1] + k[2]*l[1];
- knew[2] = k[0]*h[2] + k[1]*k[2] + k[2]*l[2];
-
- lnew[0] = l[0]*h[0] + l[1]*k[0] + l[2]*l[0];
- lnew[1] = l[0]*h[1] + l[1]*k[1] + l[2]*l[1];
- lnew[2] = l[0]*h[2] + l[1]*k[2] + l[2]*l[2];
-
- if ( is_identity(hnew, knew, lnew) ) break;
-
- memcpy(h, hnew, 3*sizeof(signed int));
- memcpy(k, knew, 3*sizeof(signed int));
- memcpy(l, lnew, 3*sizeof(signed int));
-
- }
-
- return n+1;
-}
-
-
-/* This returns the number of operations in "ops". To get the number of
- * symmetric equivalents this generates, use num_equivs() instead. */
+/* This returns the number of operations in "ops". This might be different
+ * to num_equivs() if the point group is being constructed. */
static int num_ops(const SymOpList *ops)
{
return ops->n_ops;
}
-static void update_num_equivs(SymOpList *ops)
-{
- int i, n, tot;
-
- n = num_ops(ops);
- tot = 1;
-
- for ( i=0; i<n; i++ ) {
- tot *= ops->ops[i].order;
- }
-
- ops->num_equivs = tot;
-}
-
-
/* Add a operation to a SymOpList */
static void add_symop(SymOpList *ops,
- signed int *h, signed int *k, signed int *l)
+ signed int *h, signed int *k, signed int *l,
+ int order)
{
- int n, i;
+ int n;
if ( ops->n_ops == ops->max_ops ) {
/* Pretty sure this never happens, but still... */
@@ -221,27 +190,7 @@ static void add_symop(SymOpList *ops,
ops->ops[n].h = h;
ops->ops[n].k = k;
ops->ops[n].l = l;
- ops->ops[n].order = order_of_op(h, k, l);
- ops->n_ops++;
-
- ops->divisors[0] = 1;
- for ( i=1; i<ops->n_ops; i++ ) {
- ops->divisors[i] = ops->divisors[i-1]*ops->ops[i].order;
- }
-
- update_num_equivs(ops);
-}
-
-
-static void add_copied_symop(SymOpList *ops, struct sym_op *copyme)
-{
- if ( ops->n_ops == ops->max_ops ) {
- /* Pretty sure this never happens, but still... */
- ops->max_ops += 16;
- alloc_ops(ops);
- }
-
- memcpy(&ops->ops[ops->n_ops], copyme, sizeof(*copyme));
+ ops->ops[n].order = order;
ops->n_ops++;
}
@@ -252,9 +201,9 @@ static void add_copied_symop(SymOpList *ops, struct sym_op *copyme)
* Returns: the number of equivalent reflections for a general reflection
* in point group "ops".
**/
-int num_equivs(const SymOpList *ops)
+int num_equivs(const SymOpList *ops, const SymOpMask *m)
{
- return ops->num_equivs;
+ return num_ops(ops);
}
@@ -267,22 +216,167 @@ static signed int *v(signed int h, signed int k, signed int i, signed int l)
}
+static void combine_ops(signed int *h1, signed int *k1, signed int *l1,
+ signed int *h2, signed int *k2, signed int *l2,
+ signed int *hnew, signed int *knew, signed int *lnew)
+{
+ /* Yay matrices */
+ hnew[0] = h1[0]*h2[0] + h1[1]*k2[0] + h1[2]*l2[0];
+ hnew[1] = h1[0]*h2[1] + h1[1]*k2[1] + h1[2]*l2[1];
+ hnew[2] = h1[0]*h2[2] + h1[1]*k2[2] + h1[2]*l2[2];
+
+ knew[0] = k1[0]*h2[0] + k1[1]*k2[0] + k1[2]*l2[0];
+ knew[1] = k1[0]*h2[1] + k1[1]*k2[1] + k1[2]*l2[1];
+ knew[2] = k1[0]*h2[2] + k1[1]*k2[2] + k1[2]*l2[2];
+
+ lnew[0] = l1[0]*h2[0] + l1[1]*k2[0] + l1[2]*l2[0];
+ lnew[1] = l1[0]*h2[1] + l1[1]*k2[1] + l1[2]*l2[1];
+ lnew[2] = l1[0]*h2[2] + l1[1]*k2[2] + l1[2]*l2[2];
+}
+
+
+static void expand_all_ops(signed int *hs, signed int *ks, signed int *ls,
+ int skip, SymOpList *s, SymOpList *expanded)
+{
+ int i, n;
+
+ n = num_ops(s);
+
+ for ( i=0; i<n; i++ ) {
+
+ int oi;
+ struct sym_op *opi = &s->ops[i];
+
+ for ( oi=0; oi<opi->order; oi++ ) {
+
+ int oi_it;
+ signed int *h_ordered, *k_ordered, *l_ordered;
+
+ h_ordered = malloc(3*sizeof(signed int));
+ k_ordered = malloc(3*sizeof(signed int));
+ l_ordered = malloc(3*sizeof(signed int));
+ assert(h_ordered != NULL);
+ assert(k_ordered != NULL);
+ assert(l_ordered != NULL);
+
+ memcpy(h_ordered, hs, 3*sizeof(signed int));
+ memcpy(k_ordered, ks, 3*sizeof(signed int));
+ memcpy(l_ordered, ls, 3*sizeof(signed int));
+
+ /* Apply "opi" this number of times */
+ for ( oi_it=0; oi_it<oi; oi_it++ ) {
+
+ signed int hfs[3], kfs[3], lfs[3];
+
+ combine_ops(h_ordered, k_ordered, l_ordered,
+ opi->h, opi->k, opi->l,
+ hfs, kfs, lfs);
+
+ if ( i != skip ) {
+
+ memcpy(h_ordered, hfs,
+ 3*sizeof(signed int));
+ memcpy(k_ordered, kfs,
+ 3*sizeof(signed int));
+ memcpy(l_ordered, lfs,
+ 3*sizeof(signed int));
+
+ }
+ }
+
+ STATUS("i=%i, oi=%i\n", i, oi);
+ STATUS("Creating %3i %3i %3i\n", h_ordered[0],
+ h_ordered[1],
+ h_ordered[2]);
+ STATUS(" %3i %3i %3i\n", k_ordered[0],
+ k_ordered[1],
+ k_ordered[2]);
+ STATUS(" %3i %3i %3i\n", l_ordered[0],
+ l_ordered[1],
+ l_ordered[2]);
+ STATUS("\n");
+ add_symop(expanded, h_ordered, k_ordered, l_ordered, 1);
+
+ }
+
+ }
+}
+
+
+/* Fill in the other operations for a point group starting from its
+ * generators */
+static SymOpList *expand_ops(SymOpList *s)
+{
+ int n, i;
+ SymOpList *expanded;
+
+ n = num_ops(s);
+ expanded = new_symoplist();
+ if ( expanded == NULL ) return NULL;
+ expanded->name = strdup(symmetry_name(s));
+ STATUS("%i ops to expand.\n", n);
+
+ for ( i=0; i<n; i++ ) {
+
+ int oi;
+ struct sym_op *opi = &s->ops[i];
+ STATUS("Op %i, order=%i\n", i, opi->order);
+
+ for ( oi=0; oi<opi->order; oi++ ) {
+
+ int oi_it;
+ signed int h_ordered[3], k_ordered[3], l_ordered[3];
+
+ h_ordered[0] = 1; h_ordered[1] = 0; h_ordered[2] = 0;
+ k_ordered[0] = 0; k_ordered[1] = 1; k_ordered[2] = 0;
+ l_ordered[0] = 0; l_ordered[1] = 0; l_ordered[2] = 1;
+
+ /* Apply "opi" this number of times */
+ for ( oi_it=0; oi_it<oi; oi_it++ ) {
+
+ signed int hfs[3], kfs[3], lfs[3];
+
+ combine_ops(h_ordered, k_ordered, l_ordered,
+ opi->h, opi->k, opi->l,
+ hfs, kfs, lfs);
+
+ memcpy(h_ordered, hfs, 3*sizeof(signed int));
+ memcpy(k_ordered, kfs, 3*sizeof(signed int));
+ memcpy(l_ordered, lfs, 3*sizeof(signed int));
+
+ }
+
+ STATUS("Upper: i=%i, oi=%i\n", i, oi);
+ expand_all_ops(h_ordered, k_ordered, l_ordered, i,
+ s, expanded);
+
+ }
+
+ }
+
+ free_symoplist(s);
+
+ return expanded;
+}
+
+
/********************************* 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)); /* -I */
+ add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,1), 2); /* -I */
new->name = strdup("-1");
- return new;
+ return expand_ops(new);
}
static SymOpList *make_1()
{
SymOpList *new = new_symoplist();
+ add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,1), 1); /* I */
new->name = strdup("1");
- return new;
+ return expand_ops(new);
}
@@ -291,28 +385,28 @@ static SymOpList *make_1()
static SymOpList *make_2m()
{
SymOpList *new = new_symoplist();
- add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 */
- add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* m */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 */
+ add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* m */
new->name = strdup("2/m");
- return new;
+ return expand_ops(new);
}
static SymOpList *make_2()
{
SymOpList *new = new_symoplist();
- add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 */
new->name = strdup("2");
- return new;
+ return expand_ops(new);
}
static SymOpList *make_m()
{
SymOpList *new = new_symoplist();
- add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* m */
+ add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* m */
new->name = strdup("m");
- return new;
+ return expand_ops(new);
}
@@ -321,31 +415,31 @@ static SymOpList *make_m()
static SymOpList *make_mmm()
{
SymOpList *new = new_symoplist();
- add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 */
- add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 */
- add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* m */
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 */
+ add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* m */
new->name = strdup("mmm");
- return new;
+ return expand_ops(new);
}
static SymOpList *make_222()
{
SymOpList *new = new_symoplist();
- add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 */
- add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 */
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 */
new->name = strdup("222");
- return new;
+ return expand_ops(new);
}
static SymOpList *make_mm2()
{
SymOpList *new = new_symoplist();
- add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 */
- add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* m */
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 */
+ add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* m */
new->name = strdup("mm2");
- return new;
+ return expand_ops(new);
}
@@ -360,7 +454,7 @@ static SymOpList *make_4m()
static SymOpList *make_4()
{
SymOpList *new = new_symoplist();
- add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1)); /* 4 */
+ add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 */
return NULL;
}
@@ -488,8 +582,8 @@ static void do_op(const struct sym_op *op,
signed int *he, signed int *ke, signed int *le)
{
*he = h*op->h[0] + k*op->h[1] + l*op->h[2];
- *ke = h*op->k[0] + k*op->h[1] + l*op->k[2];
- *le = h*op->l[0] + k*op->h[1] + l*op->l[2];
+ *ke = h*op->k[0] + k*op->k[1] + l*op->k[2];
+ *le = h*op->l[0] + k*op->l[1] + l*op->l[2];
}
@@ -513,37 +607,31 @@ static void do_op(const struct sym_op *op,
* given reflection is a special high-symmetry one), call special_position()
* first to get a "specialised" SymOpList and use that instead.
**/
-void get_equiv(const SymOpList *ops, int idx,
+void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx,
signed int h, signed int k, signed int l,
signed int *he, signed int *ke, signed int *le)
{
- int sig[32];
- int i, n, r;
+ int i, n;
n = num_ops(ops);
- r = idx;
- for ( i=n-1; i>=0; i-- ) {
- sig[i] = r / ops->divisors[i];
- r = r % ops->divisors[i];
- assert(sig[i] < ops->ops[i].order);
+ for ( i=idx; i<n; i++ ) {
+ if ( (m == NULL) || m->mask[i] ) {
+ do_op(&ops->ops[i], h, k, l, he, ke, le);
+ return;
+ }
}
- for ( i=0; i<n; i++ ) {
-
- int s;
+ ERROR("Index %i out of range for point group '%s'\n", idx,
+ symmetry_name(ops));
- /* Do this operation "sig[i]" times */
- for ( s=0; s<sig[i]; s++ ) {
- do_op(&ops->ops[i], h, k, l, &h, &k, &l);
- }
-
- }
+ *he = 0; *ke = 0; *le = 0;
}
/**
* special_position:
* @ops: A %SymOpList, usually corresponding to a point group
+ * @m: A %SymOpMask created with new_symopmask()
* @h: index of a reflection
* @k: index of a reflection
* @l: index of a reflection
@@ -554,26 +642,47 @@ void get_equiv(const SymOpList *ops, int idx,
*
* Returns: the "specialised" %SymOpList.
**/
-SymOpList *special_position(const SymOpList *ops,
- signed int h, signed int k, signed int l)
+void special_position(const SymOpList *ops, SymOpMask *m,
+ signed int h, signed int k, signed int l)
{
int i, n;
- SymOpList *specialised;
+ signed int *htest;
+ signed int *ktest;
+ signed int *ltest;
- n = num_ops(ops);
- specialised = new_symoplist();
+ assert(m->list = ops);
+
+ n = num_equivs(ops, NULL);
+ htest = malloc(n*sizeof(signed int));
+ ktest = malloc(n*sizeof(signed int));
+ ltest = malloc(n*sizeof(signed int));
for ( i=0; i<n; i++ ) {
- signed int ht, kt, lt;
+ signed int he, ke, le;
+ int j;
+
+ get_equiv(ops, NULL, i, h, k, l, &he, &ke, &le);
+
+ m->mask[i] = 1;
+ for ( j=0; j<i; j++ ) {
+ if ( (he==htest[j]) && (ke==ktest[j])
+ && (le==ltest[j]) )
+ {
+ m->mask[i] = 0;
+ break; /* Only need to find one */
+ }
+ }
- do_op(&ops->ops[i], h, k, l, &ht, &kt, &lt);
- if ( (h==ht) || (k==kt) || (l==lt) ) continue;
- add_copied_symop(specialised, &ops->ops[i]);
+ htest[i] = he;
+ ktest[i] = ke;
+ ltest[i] = le;
}
- return specialised;
+ free(htest);
+ free(ktest);
+ free(ltest);
}
@@ -585,12 +694,12 @@ void get_asymm(const SymOpList *ops,
int p;
signed int best_h, best_k, best_l;
- nequiv = num_equivs(ops);
+ nequiv = num_equivs(ops, NULL);
best_h = h; best_k = k; best_l = l;
for ( p=0; p<nequiv; p++ ) {
- get_equiv(ops, p, h, k, l, hp, kp, lp);
+ get_equiv(ops, NULL, p, h, k, l, hp, kp, lp);
if ( h > best_h ) {
best_h = h; best_k = k; best_l = l;
diff --git a/src/symmetry.h b/src/symmetry.h
index 4dcaae06..45ffe318 100644
--- a/src/symmetry.h
+++ b/src/symmetry.h
@@ -24,22 +24,27 @@
**/
typedef struct _symoplist SymOpList;
-extern void free_symoplist(SymOpList *ops);
+/**
+ * SymOpMask
+ *
+ * Opaque type.
+ **/
+typedef struct _symopmask SymOpMask;
+extern void free_symoplist(SymOpList *ops);
extern SymOpList *get_pointgroup(const char *sym);
-
extern const char *symmetry_name(const SymOpList *ops);
-extern SymOpList *special_position(const SymOpList *ops,
- signed int h, signed int k, signed int l);
+extern SymOpMask *new_symopmask(const SymOpList *ops);
+extern void free_symopmask(SymOpMask *m);
+extern void special_position(const SymOpList *ops, SymOpMask *m,
+ signed int h, signed int k, signed int l);
extern void get_asymm(const SymOpList *ops,
signed int h, signed int k, signed int l,
signed int *hp, signed int *kp, signed int *lp);
-
-extern int num_equivs(const SymOpList *ops);
-
-extern void get_equiv(const SymOpList *ops, int idx,
+extern int num_equivs(const SymOpList *ops, const SymOpMask *m);
+extern void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx,
signed int h, signed int k, signed int l,
signed int *he, signed int *ke, signed int *le);