diff options
Diffstat (limited to 'src/symmetry.c')
| -rw-r--r-- | src/symmetry.c | 464 |
1 files changed, 464 insertions, 0 deletions
diff --git a/src/symmetry.c b/src/symmetry.c new file mode 100644 index 0000000..65aef98 --- /dev/null +++ b/src/symmetry.c @@ -0,0 +1,464 @@ +/* + * symmetry.c + * + * Symmetry stuff + * + * (c) 2006-2008 Thomas White <taw27@cam.ac.uk> + * + * synth2d - Two-Dimensional Crystallographic Fourier Synthesis + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <math.h> +#include <assert.h> + +#include "symmetry.h" +#include "reflist.h" +#include "model.h" + +#define is_odd(a) ((a)%2==1) + +/* Return a list of atoms containing the given atom and all its symmetry equivalents */ +AtomicModel *symmetry_generate_equivalent_atoms(AtomicModel *model, size_t j, ModelTarget target) { + + AtomicModel *list; + + list = model_new(); + list->atoms[list->n_atoms] = model->atoms[j]; list->n_atoms++; + + if ( target == MODEL_TARGET_DISPLAY ) { + + if ( model->atoms[j].x == 0 ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = 1.0; + list->n_atoms++; + } + + if ( model->atoms[j].y == 0 ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].y = 1.0; + list->n_atoms++; + } + + if ( (model->atoms[j].x == 0) && (model->atoms[j].y == 0) ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = 1.0; + list->atoms[list->n_atoms].y = 1.0; + list->n_atoms++; + } + + } + + if ( model->sym & SYMMETRY_CENTRE ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = 1-list->atoms[list->n_atoms].x; + list->atoms[list->n_atoms].y = 1-list->atoms[list->n_atoms].y; + list->n_atoms++; + } + + if ( model->sym & SYMMETRY_MIRROR_HORIZONTAL ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].y = 1-list->atoms[list->n_atoms].y; + list->n_atoms++; + } + + if ( model->sym & SYMMETRY_MIRROR_VERTICAL ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = 1-list->atoms[list->n_atoms].x; + list->n_atoms++; + } + + if ( model->sym & SYMMETRY_GLIDE_HORIZONTAL ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = fmod(0.5 + list->atoms[list->n_atoms].x, 1); + list->atoms[list->n_atoms].y = 1-list->atoms[list->n_atoms].y; + list->n_atoms++; + } + + if ( model->sym & SYMMETRY_GLIDE_VERTICAL ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = 1-list->atoms[list->n_atoms].x; + list->atoms[list->n_atoms].y = fmod(0.5 + list->atoms[list->n_atoms].y, 1); + list->n_atoms++; + } + + if ( model->sym & SYMMETRY_GLIDE_HORIZONTAL_QUARTER ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = fmod(0.5 + list->atoms[list->n_atoms].x, 1); + list->atoms[list->n_atoms].y = fmod(1.5 - list->atoms[list->n_atoms].y, 1); + list->n_atoms++; + } + + if ( model->sym & SYMMETRY_GLIDE_VERTICAL_QUARTER ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = fmod(1.5 - list->atoms[list->n_atoms].x, 1); + list->atoms[list->n_atoms].y = fmod(0.5 + list->atoms[list->n_atoms].y, 1); + list->n_atoms++; + } + + /* Fix cm */ + if ( (model->sym & SYMMETRY_GLIDE_HORIZONTAL_QUARTER) && (model->sym & SYMMETRY_MIRROR_HORIZONTAL) ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = fmod(0.5 + list->atoms[list->n_atoms].x, 1); + list->atoms[list->n_atoms].y = 1 - fmod(1.5 - list->atoms[list->n_atoms].y, 1); + list->n_atoms++; + } else if ( (model->sym & SYMMETRY_GLIDE_VERTICAL_QUARTER) && (model->sym & SYMMETRY_MIRROR_VERTICAL) ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = 1 - fmod(1.5 - list->atoms[list->n_atoms].x, 1); + list->atoms[list->n_atoms].y = fmod(0.5 + list->atoms[list->n_atoms].y, 1); + list->n_atoms++; + } + + /* Fix c2mm */ + if ( (model->sym & SYMMETRY_GLIDE_VERTICAL_QUARTER) && (model->sym & SYMMETRY_MIRROR_HORIZONTAL) ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = fmod(1.5 - list->atoms[list->n_atoms].x, 1); + list->atoms[list->n_atoms].y = 1 - fmod(0.5 + list->atoms[list->n_atoms].y, 1); + list->n_atoms++; + } + + if ( model->sym & SYMMETRY_MIRROR_HORIZONTAL_QUARTER ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = list->atoms[list->n_atoms].x; + list->atoms[list->n_atoms].y = fmod(1.5 - list->atoms[list->n_atoms].y, 1); + list->n_atoms++; + } + + if ( model->sym & SYMMETRY_MIRROR_VERTICAL_QUARTER ) { + list->atoms[list->n_atoms] = model->atoms[j]; + list->atoms[list->n_atoms].x = fmod(1.5 - list->atoms[list->n_atoms].x, 1); + list->atoms[list->n_atoms].y = list->atoms[list->n_atoms].y; + list->n_atoms++; + } + + return list; + +} + +ReflectionList *symmetry_generate_equivalent_reflections(Symmetry sym, signed int h, signed int k, signed int l) { + + ReflectionList *result; + + result = reflist_new(); + + /* Remove the minus in front of 'l' to make even HOLZ conditional transforms real. + * This is generally the wrong thing to do. */ + if ( sym & SYMMETRY_FRIEDEL ) reflist_addref_deltatheta(result, -h, -k, -l, 0, -1); + + if ( sym & SYMMETRY_MIRROR_HORIZONTAL ) reflist_addref_deltatheta(result, h, -k, l, 0, 1); + if ( sym & SYMMETRY_MIRROR_VERTICAL ) reflist_addref_deltatheta(result, -h, k, l, 0, 1); + if ( sym & SYMMETRY_MIRROR_DIAGONAL ) reflist_addref_deltatheta(result, -h, -k, l, 0, 1); + + if ( sym & SYMMETRY_GLIDE_HORIZONTAL ) reflist_addref_deltatheta(result, h, -k, l, M_PI, 1); + if ( sym & SYMMETRY_GLIDE_VERTICAL ) reflist_addref_deltatheta(result, -h, k, l, M_PI, 1); + + if ( sym & SYMMETRY_GLIDE_HORIZONTAL_QUARTER ) { + if ( is_odd(abs(h)+abs(k)) ) { + reflist_addref_deltatheta(result, h, -k, l, M_PI, 1); + } else { + reflist_addref_deltatheta(result, h, -k, l, 0, 1); + } + } + + if ( sym & SYMMETRY_GLIDE_VERTICAL_QUARTER ) { + if ( is_odd(abs(h)+abs(k)) ) { + reflist_addref_deltatheta(result, -h, k, l, M_PI, 1); + } else { + reflist_addref_deltatheta(result, -h, k, l, 0, 1); + } + } + + if ( sym & SYMMETRY_ROTATION_4 ) { + if ( sym & SYMMETRY_GLIDE_DIAGONAL ) { + if ( is_odd(abs(h)+abs(k)) ) { + reflist_addref_deltatheta(result, h, -k, l, M_PI, 1); + reflist_addref_deltatheta(result, -h, k, l, M_PI, 1); + reflist_addref_deltatheta(result, -h, -k, l, 0, 1); + reflist_addref_deltatheta(result, k, h, l, 0, 1); + reflist_addref_deltatheta(result, -k, h, l, M_PI, 1); + reflist_addref_deltatheta(result, k, -h, l, M_PI, 1); + reflist_addref_deltatheta(result, -k, -h, l, 0, 1); + } else { + reflist_addref_deltatheta(result, h, -k, l, 0, 1); + reflist_addref_deltatheta(result, -h, k, l, 0, 1); + reflist_addref_deltatheta(result, -h, -k, l, 0, 1); + reflist_addref_deltatheta(result, k, h, l, 0, 1); + reflist_addref_deltatheta(result, -k, h, l, 0, 1); + reflist_addref_deltatheta(result, k, -h, l, 0, 1); + reflist_addref_deltatheta(result, -k, -h, l, 0, 1); + } + } else { + reflist_addref_deltatheta(result, -h, -k, l, 0, 1); + reflist_addref_deltatheta(result, -h, k, l, 0, 1); + reflist_addref_deltatheta(result, h, -k, l, 0, 1); + } + } + + return result; + +} + +static double symmetry_realphase(double phk) { + + if ( (phk >= 0) && (phk < M_PI_2) ) phk = 0; + else if ( (phk > M_PI_2) && (phk < M_PI) ) phk = M_PI; + else if ( (phk < 0) && (phk > -M_PI/2) ) phk = 0; + else if ( (phk <= -M_PI_2) && (phk > -M_PI) ) phk = M_PI; + + return phk; + +} + +void symmetry_centricity(ReflectionList *reflections, unsigned int i, Symmetry sym, SymFlags flags) { + + double am, ph; + signed int h, k; + + ph = 0; + h = reflections->refs[i].h; + k = reflections->refs[i].k; + + am = reflections->refs[i].amplitude; + if ( flags & SYMFLAG_PHASES_KNOWN ) ph = fmod(reflections->refs[i].phase_known, 2*M_PI); + if ( flags & SYMFLAG_PHASES_CALC ) ph = fmod(reflections->refs[i].phase_calc, 2*M_PI); + + if ( (sym & SYMMETRY_MIRROR_HORIZONTAL) && (sym & SYMMETRY_FRIEDEL) && (h == 0) ) { + ph = symmetry_realphase(ph); + } + + if ( (sym & SYMMETRY_MIRROR_VERTICAL) && (sym & SYMMETRY_FRIEDEL) && (k == 0) ) { + ph = symmetry_realphase(ph); + } + + if ( (sym & SYMMETRY_CENTRE) && (sym & SYMMETRY_FRIEDEL) ) { + ph = symmetry_realphase(ph); + } + + if (( ( (sym & SYMMETRY_GLIDE_HORIZONTAL) && (sym & SYMMETRY_MIRROR_VERTICAL) ) + || ( (sym & SYMMETRY_GLIDE_HORIZONTAL_QUARTER) && (sym & SYMMETRY_MIRROR_VERTICAL) ) + || ( (sym & SYMMETRY_GLIDE_VERTICAL) && (sym & SYMMETRY_MIRROR_HORIZONTAL) ) + || ( (sym & SYMMETRY_GLIDE_VERTICAL_QUARTER) && (sym & SYMMETRY_MIRROR_HORIZONTAL) ) ) + && (is_odd(abs(h)+abs(k))) ) { + if ( (ph >= 0) && (ph < M_PI) ) ph = M_PI_2; + else if ( (ph >= M_PI_2) && (ph < 2*M_PI) ) ph = -M_PI_2; + else if ( (ph < 0) && (ph > -M_PI) ) ph = -M_PI_2; + else if ( (ph <= -M_PI) ) ph = M_PI_2; + } + + if ( flags & SYMFLAG_PHASES_KNOWN ) reflections->refs[i].phase_known = ph; + if ( flags & SYMFLAG_PHASES_CALC ) reflections->refs[i].phase_calc = ph; + +} + +unsigned int symmetry_reflection_allowed(signed int h, signed int k, signed int l, Symmetry sym) { + + if ( (sym & SYMMETRY_GLIDE_HORIZONTAL) && (h == 0) && (k % 2) ) return 0; + if ( (sym & SYMMETRY_GLIDE_HORIZONTAL_QUARTER) && (h == 0) && (k % 2) ) return 0; + if ( (sym & SYMMETRY_GLIDE_VERTICAL) && (k == 0) && (h % 2) ) return 0; + if ( (sym & SYMMETRY_GLIDE_VERTICAL_QUARTER) && (k == 0) && (h % 2) ) return 0; + if ( (sym & SYMMETRY_ROTATION_4 ) && ( sym & SYMMETRY_GLIDE_DIAGONAL ) && (k == 0) && (h % 2) ) return 0; + if ( (sym & SYMMETRY_ROTATION_4 ) && ( sym & SYMMETRY_GLIDE_DIAGONAL ) && (h == 0) && (k % 2) ) return 0; + + return 1; + +} + +/* Symmetrise a list of reflections (expanding all equivalents to P1) */ +double symmetry_symmetrise(ReflectionList *reflections, Symmetry sym, SymFlags flags) { + + unsigned int i = 0; + double r_sym = 0; + double r_sym_norm = 0; + unsigned int new_ref = 0; + unsigned int elim_ref = 0; + double am_elim = 0; + unsigned int n = reflections->n_reflections; + + for ( i=1; i<n; i++ ) { /* Don't symmetrise 000 */ + + signed int h, k, l; + + h = reflections->refs[i].h; + k = reflections->refs[i].k; + l = reflections->refs[i].l; + + if ( (h==0) && (k==0) && (l==0) ) continue; + if ( symmetry_reflection_allowed(h, k, l, sym) ) { + + double ph = 69; + unsigned int j; + double av; + unsigned int av_n; + ReflectionList *equivalents; + + equivalents = symmetry_generate_equivalent_reflections(sym, h, k, l); + symmetry_centricity(reflections, i, sym, flags); + + /* First pass: calculate average amplitude, determine phase */ + av = reflections->refs[i].amplitude; av_n = 1; + for ( j=1; j<equivalents->n_reflections; j++ ) { + unsigned int p; + p = reflist_inlist(reflections, equivalents->refs[j].h, equivalents->refs[j].k, equivalents->refs[j].l); + if ( p ) { + av += reflections->refs[p].amplitude; + av_n++; + } + } + av = av / av_n; + r_sym += fabs(av - reflections->refs[i].amplitude); + r_sym_norm += reflections->refs[i].amplitude; + if ( flags & SYMFLAG_PHASES_KNOWN ) ph = fmod(reflections->refs[i].phase_known, 2*M_PI); + if ( flags & SYMFLAG_PHASES_CALC ) ph = fmod(reflections->refs[i].phase_calc, 2*M_PI); + //printf("Phase of %3i %3i %3i is %f\n", h, k, l, ph); + + /* Second pass: set all equivalent reflections to the average amplitude */ + reflections->refs[i].amplitude = av; + for ( j=1; j<equivalents->n_reflections; j++ ) { + + unsigned int p; + double phn; + p = reflist_inlist(reflections, equivalents->refs[j].h, equivalents->refs[j].k, equivalents->refs[j].l); + if ( p ) { + r_sym += fabs(av - reflections->refs[p].amplitude); + r_sym_norm += reflections->refs[p].amplitude; + reflections->refs[p].amplitude = av; + } else { + p = reflist_addref_am(reflections, equivalents->refs[j].h, equivalents->refs[j].k, equivalents->refs[j].l, av); + //printf("SY: Generating %3i %3i %3i am=%f\n", equivalents->refs[j].h, + // equivalents->refs[j].k, equivalents->refs[j].l, av); + new_ref++; + } + phn = (ph + equivalents->refs[j].delta_theta) * equivalents->refs[j].multiplier; + if ( flags & SYMFLAG_PHASES_KNOWN ) { + reflections->refs[p].phase_known = phn; + reflections->refs[p].phase_known_set = 1; + } + if ( flags & SYMFLAG_PHASES_CALC ) { + reflections->refs[p].phase_calc = phn; + reflections->refs[p].phase_calc_set = 1; + } + //printf("Set phase of %3i %3i %3i to %f (dt=%f, mul=%i)\n", equivalents->refs[j].h, + //equivalents->refs[j].k, equivalents->refs[j].l, + //phn, equivalents->refs[j].delta_theta, equivalents->refs[j].multiplier); + + + } + + reflist_free(equivalents); + + } else { + + am_elim += reflections->refs[i].amplitude; + reflist_delref(reflections, h, k, l); + //printf("SY: Eliminating systematically absent reflection %i %i %i\n", h, k, l); + elim_ref++; + i--; + + } + + } + + if ( r_sym > 0 ) { + r_sym = r_sym / r_sym_norm; + } + //printf("SY: R_sym = %.2f%%, %i reflections generated, %i reflections eliminated (total amplitude %f)\n", r_sym*100, new_ref, elim_ref, am_elim); + + return r_sym; + +} + +Symmetry symmetry_encode(const char *symmetry) { + + if ( strcmp(symmetry, "p1") == 0 ) return PLANEGROUP_P1; + if ( strcmp(symmetry, "p2") == 0 ) return PLANEGROUP_P2; + if ( strcmp(symmetry, "pm (m // x)") == 0 ) return PLANEGROUP_PM_X; + if ( strcmp(symmetry, "pm (m // y)") == 0 ) return PLANEGROUP_PM_Y; + if ( strcmp(symmetry, "pg (g // x)") == 0 ) return PLANEGROUP_PG_X; + if ( strcmp(symmetry, "pg (g // y)") == 0 ) return PLANEGROUP_PG_Y; + if ( strcmp(symmetry, "cm (m // x)") == 0 ) return PLANEGROUP_CM_X; + if ( strcmp(symmetry, "cm (m // y)") == 0 ) return PLANEGROUP_CM_Y; + if ( strcmp(symmetry, "p2mm") == 0 ) return PLANEGROUP_P2MM; + if ( strcmp(symmetry, "p2mg (m // x)") == 0 ) return PLANEGROUP_P2MG_X; + if ( strcmp(symmetry, "p2mg (m // y)") == 0 ) return PLANEGROUP_P2MG_Y; + if ( strcmp(symmetry, "p2gg") == 0 ) return PLANEGROUP_P2GG; + if ( strcmp(symmetry, "c2mm") == 0 ) return PLANEGROUP_C2MM; + if ( strcmp(symmetry, "p4") == 0 ) return PLANEGROUP_P4; + if ( strcmp(symmetry, "p4mm") == 0 ) return PLANEGROUP_P4MM; + if ( strcmp(symmetry, "p4gm") == 0 ) return PLANEGROUP_P4GM; + if ( strcmp(symmetry, "p3") == 0 ) return PLANEGROUP_P3; + if ( strcmp(symmetry, "p3m1") == 0 ) return PLANEGROUP_P3M1; + if ( strcmp(symmetry, "p31m") == 0 ) return PLANEGROUP_P31M; + if ( strcmp(symmetry, "p6") == 0 ) return PLANEGROUP_P6; + if ( strcmp(symmetry, "p6mm") == 0 ) return PLANEGROUP_P6MM; + + fprintf(stderr, "Unrecognised symmetry identifier '%s'\n", symmetry); + return SYMMETRY_IDENTITY; + +} + +const char *symmetry_decode(Symmetry sym) { + + if ( sym == PLANEGROUP_P1 ) return "p1"; + if ( sym == PLANEGROUP_P2 ) return "p2"; + if ( sym == PLANEGROUP_PM_X ) return "pm (m // x)"; + if ( sym == PLANEGROUP_PM_Y ) return "pm (m // y)"; + if ( sym == PLANEGROUP_PG_X ) return "pg (g // x)"; + if ( sym == PLANEGROUP_PG_Y ) return "pg (g // y)"; + if ( sym == PLANEGROUP_CM_X ) return "cm (m // x)"; + if ( sym == PLANEGROUP_CM_Y ) return "cm (m // y)"; + if ( sym == PLANEGROUP_P2MM ) return "p2mm"; + if ( sym == PLANEGROUP_P2MG_X ) return "p2mg (m // x)"; + if ( sym == PLANEGROUP_P2MG_Y ) return "p2mg (m // y)"; + if ( sym == PLANEGROUP_P2GG ) return "p2gg"; + if ( sym == PLANEGROUP_C2MM ) return "c2mm"; + if ( sym == PLANEGROUP_P4 ) return "p4"; + if ( sym == PLANEGROUP_P4MM ) return "p4mm"; + if ( sym == PLANEGROUP_P4GM ) return "p4gm"; + if ( sym == PLANEGROUP_P3 ) return "p3"; + if ( sym == PLANEGROUP_P3M1 ) return "p3m1"; + if ( sym == PLANEGROUP_P31M ) return "p31m"; + if ( sym == PLANEGROUP_P6 ) return "p6"; + if ( sym == PLANEGROUP_P6MM ) return "p6mm"; + if ( sym == (PLANEGROUP_P1 | SYMMETRY_FRIEDEL) ) return "p1 & Friedel"; + if ( sym == (PLANEGROUP_P2 | SYMMETRY_FRIEDEL) ) return "p2 & Friedel"; + if ( sym == (PLANEGROUP_PM_X | SYMMETRY_FRIEDEL) ) return "pm (m // x) & Friedel"; + if ( sym == (PLANEGROUP_PM_Y | SYMMETRY_FRIEDEL) ) return "pm (m // y) & Friedel"; + if ( sym == (PLANEGROUP_PG_X | SYMMETRY_FRIEDEL) ) return "pg (g // x) & Friedel"; + if ( sym == (PLANEGROUP_PG_Y | SYMMETRY_FRIEDEL) ) return "pg (g // y) & Friedel"; + if ( sym == (PLANEGROUP_CM_X | SYMMETRY_FRIEDEL) ) return "cm (m // x) & Friedel"; + if ( sym == (PLANEGROUP_CM_Y | SYMMETRY_FRIEDEL) ) return "cm (m // y) & Friedel"; + if ( sym == (PLANEGROUP_P2MM | SYMMETRY_FRIEDEL) ) return "p2mm & Friedel"; + if ( sym == (PLANEGROUP_P2MG_X | SYMMETRY_FRIEDEL) ) return "p2mg (m // x) & Friedel"; + if ( sym == (PLANEGROUP_P2MG_Y | SYMMETRY_FRIEDEL) ) return "p2mg (m // y) & Friedel"; + if ( sym == (PLANEGROUP_P2GG | SYMMETRY_FRIEDEL) ) return "p2gg & Friedel"; + if ( sym == (PLANEGROUP_C2MM | SYMMETRY_FRIEDEL) ) return "c2mm & Friedel"; + if ( sym == (PLANEGROUP_P4 | SYMMETRY_FRIEDEL) ) return "p4 & Friedel"; + if ( sym == (PLANEGROUP_P4MM | SYMMETRY_FRIEDEL) ) return "p4mm & Friedel"; + if ( sym == (PLANEGROUP_P4GM | SYMMETRY_FRIEDEL) ) return "p4gm & Friedel"; + if ( sym == (PLANEGROUP_P3 | SYMMETRY_FRIEDEL) ) return "p3 & Friedel"; + if ( sym == (PLANEGROUP_P3M1 | SYMMETRY_FRIEDEL) ) return "p3m1 & Friedel"; + if ( sym == (PLANEGROUP_P31M | SYMMETRY_FRIEDEL) ) return "p31m & Friedel"; + if ( sym == (PLANEGROUP_P6 | SYMMETRY_FRIEDEL) ) return "p6 & Friedel"; + if ( sym == (PLANEGROUP_P6MM | SYMMETRY_FRIEDEL) ) return "p6mm & Friedel"; + + return "(unknown symmetry)"; + +} + +void symmetry_symmetrise_array(fftw_complex *in, signed int width, signed int height, Symmetry sym) { + + ReflectionList *reflections; + + reflections = reflist_new_from_array(in, width, height); + symmetry_symmetrise(reflections, sym, SYMFLAG_PHASES_KNOWN); + reflist_fill_array(in, reflections, width, height); + reflist_free(reflections); + +} |