From b143429764665a75dd3baf8c5115bf8553d18d71 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 20 Jul 2011 17:53:27 +0200 Subject: Symmetry stuff --- Makefile.am | 9 +- Makefile.in | 28 +++- src/diffraction.c | 8 +- src/get_hkl.c | 29 +++- src/povray.c | 22 ++- src/powder_plot.c | 10 +- src/process_hkl.c | 3 +- src/reflist-utils.c | 23 ++- src/render_hkl.c | 14 +- src/symmetry.c | 423 +++++++++++++++++++++++++++++++++------------------- src/symmetry.h | 21 ++- 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(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 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 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; ilist = list; + m->mask = malloc(sizeof(int)*list->n_ops); + if ( m->mask == NULL ) { + free(m); + return NULL; + } + + for ( i=0; in_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; iops[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; in_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; iops[i]; + + for ( oi=0; oiorder; 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_ith, 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; iops[i]; + STATUS("Op %i, order=%i\n", i, opi->order); + + for ( oi=0; oiorder; 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_ith, 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; imask[i] ) { + do_op(&ops->ops[i], h, k, l, he, ke, le); + return; + } } - for ( i=0; iops[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; imask[i] = 1; + for ( j=0; jmask[i] = 0; + break; /* Only need to find one */ + } + } - do_op(&ops->ops[i], h, k, l, &ht, &kt, <); - 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 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); -- cgit v1.2.3