diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/Makefile.am | 2 | ||||
-rw-r--r-- | src/Makefile.in | 7 | ||||
-rw-r--r-- | src/diffraction-gpu.c | 3 | ||||
-rw-r--r-- | src/diffraction-gpu.h | 6 | ||||
-rw-r--r-- | src/diffraction.c | 122 | ||||
-rw-r--r-- | src/diffraction.h | 5 | ||||
-rw-r--r-- | src/indexamajig.c | 52 | ||||
-rw-r--r-- | src/pattern_sim.c | 31 | ||||
-rw-r--r-- | src/utils.h | 5 |
9 files changed, 185 insertions, 48 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index 6031e441..98508dd0 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -11,7 +11,7 @@ AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c cell.c \ hdf5-file.c detector.c sfac.c peaks.c reflections.c \ - beam-parameters.c + beam-parameters.c symmetry.c if HAVE_OPENCL pattern_sim_SOURCES += diffraction-gpu.c cl-utils.c endif diff --git a/src/Makefile.in b/src/Makefile.in index a427a4ef..4cc04b36 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -119,12 +119,13 @@ indexamajig_OBJECTS = $(am_indexamajig_OBJECTS) indexamajig_DEPENDENCIES = am__pattern_sim_SOURCES_DIST = pattern_sim.c diffraction.c utils.c \ image.c cell.c hdf5-file.c detector.c sfac.c peaks.c \ - reflections.c beam-parameters.c diffraction-gpu.c cl-utils.c + reflections.c beam-parameters.c symmetry.c diffraction-gpu.c \ + cl-utils.c am_pattern_sim_OBJECTS = pattern_sim.$(OBJEXT) diffraction.$(OBJEXT) \ utils.$(OBJEXT) image.$(OBJEXT) cell.$(OBJEXT) \ hdf5-file.$(OBJEXT) detector.$(OBJEXT) sfac.$(OBJEXT) \ peaks.$(OBJEXT) reflections.$(OBJEXT) \ - beam-parameters.$(OBJEXT) $(am__objects_1) + beam-parameters.$(OBJEXT) symmetry.$(OBJEXT) $(am__objects_1) pattern_sim_OBJECTS = $(am_pattern_sim_OBJECTS) pattern_sim_DEPENDENCIES = am_powder_plot_OBJECTS = powder_plot.$(OBJEXT) cell.$(OBJEXT) \ @@ -284,7 +285,7 @@ AM_CFLAGS = -Wall AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c \ cell.c hdf5-file.c detector.c sfac.c peaks.c reflections.c \ - beam-parameters.c $(am__append_2) + beam-parameters.c symmetry.c $(am__append_2) pattern_sim_LDADD = @LIBS@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \ reflections.c symmetry.c stream.c beam-parameters.c diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c index 324cc4cd..2d8eda95 100644 --- a/src/diffraction-gpu.c +++ b/src/diffraction-gpu.c @@ -330,7 +330,8 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, /* Setup the OpenCL stuff, create buffers, load the structure factor table */ struct gpu_context *setup_gpu(int no_sfac, struct image *image, - const double *intensities, int dev_num) + const double *intensities, unsigned char *flags, + int dev_num) { struct gpu_context *gctx; cl_uint nplat; diff --git a/src/diffraction-gpu.h b/src/diffraction-gpu.h index 4a5815bf..7e6048b3 100644 --- a/src/diffraction-gpu.h +++ b/src/diffraction-gpu.h @@ -26,7 +26,8 @@ struct gpu_context; extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, int na, int nb, int nc, UnitCell *ucell); extern struct gpu_context *setup_gpu(int no_sfac, struct image *image, - const double *intensities, int dev_num); + const double *intensities, + const unsigned char *flags, int dev_num); extern void cleanup_gpu(struct gpu_context *gctx); #else @@ -39,7 +40,8 @@ static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, } static struct gpu_context *setup_gpu(int no_sfac, struct image *image, - const double *intensities, int dev_num) + const double *intensities, + const unsigned char *flags, int dev_num) { return NULL; } diff --git a/src/diffraction.c b/src/diffraction.c index 93ac7f8e..aa06c44b 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -23,6 +23,7 @@ #include "diffraction.h" #include "sfac.h" #include "beam-parameters.h" +#include "symmetry.h" #define SAMPLING (4) @@ -69,8 +70,63 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az, } -static double interpolate_linear(const double *ref, - float hd, signed int k, signed int l) +static double sym_lookup_intensity(const double *intensities, + const unsigned char *flags, const char *sym, + signed int h, signed int k, signed int l) +{ + int i; + double ret = 0.0; + + for ( i=0; i<num_general_equivs(sym); i++ ) { + + signed int he; + signed int ke; + signed int le; + double f, val; + + get_general_equiv(h, k, l, &he, &ke, &le, sym, i); + + f = (double)lookup_flag(flags, he, ke, le); + val = lookup_intensity(intensities, he, ke, le); + + ret += f*val; + + } + + return ret; +} + + +static double sym_lookup_phase(const double *phases, + const unsigned char *flags, const char *sym, + signed int h, signed int k, signed int l) +{ + int i; + double ret = 0.0; + + for ( i=0; i<num_general_equivs(sym); i++ ) { + + signed int he; + signed int ke; + signed int le; + double f, val; + + get_general_equiv(h, k, l, &he, &ke, &le, sym, i); + + f = (double)lookup_flag(flags, he, ke, le); + val = lookup_phase(phases, he, ke, le); + + ret += f*val; + + } + + return ret; +} + + +static double interpolate_linear(const double *ref, const unsigned char *flags, + const char *sym, float hd, + signed int k, signed int l) { signed int h; double val1, val2; @@ -81,8 +137,8 @@ static double interpolate_linear(const double *ref, f = hd - (float)h; assert(f >= 0.0); - val1 = lookup_intensity(ref, h, k, l); - val2 = lookup_intensity(ref, h+1, k, l); + val1 = sym_lookup_intensity(ref, flags, sym, h, k, l); + val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l); val1 = val1; val2 = val2; @@ -92,6 +148,7 @@ static double interpolate_linear(const double *ref, static double interpolate_bilinear(const double *ref, + const unsigned char *flags, const char *sym, float hd, float kd, signed int l) { signed int k; @@ -103,14 +160,15 @@ static double interpolate_bilinear(const double *ref, f = kd - (float)k; assert(f >= 0.0); - val1 = interpolate_linear(ref, hd, k, l); - val2 = interpolate_linear(ref, hd, k+1, l); + val1 = interpolate_linear(ref, flags, sym, hd, k, l); + val2 = interpolate_linear(ref, flags, sym, hd, k+1, l); return (1.0-f)*val1 + f*val2; } static double interpolate_intensity(const double *ref, + const unsigned char *flags, const char *sym, float hd, float kd, float ld) { signed int l; @@ -122,8 +180,8 @@ static double interpolate_intensity(const double *ref, f = ld - (float)l; assert(f >= 0.0); - val1 = interpolate_bilinear(ref, hd, kd, l); - val2 = interpolate_bilinear(ref, hd, kd, l+1); + val1 = interpolate_bilinear(ref, flags, sym, hd, kd, l); + val2 = interpolate_bilinear(ref, flags, sym, hd, kd, l+1); return (1.0-f)*val1 + f*val2; } @@ -131,6 +189,8 @@ static double interpolate_intensity(const double *ref, static double complex interpolate_phased_linear(const double *ref, const double *phases, + const unsigned char *flags, + const char *sym, float hd, signed int k, signed int l) { @@ -146,10 +206,10 @@ static double complex interpolate_phased_linear(const double *ref, f = hd - (float)h; assert(f >= 0.0); - val1 = lookup_intensity(ref, h, k, l); - val2 = lookup_intensity(ref, h+1, k, l); - ph1 = lookup_phase(phases, h, k, l); - ph2 = lookup_phase(phases, h+1, k, l); + val1 = sym_lookup_intensity(ref, flags, sym, h, k, l); + val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l); + ph1 = sym_lookup_phase(phases, flags, sym, h, k, l); + ph2 = sym_lookup_phase(phases, flags, sym, h+1, k, l); val1 = val1; val2 = val2; @@ -169,6 +229,8 @@ static double complex interpolate_phased_linear(const double *ref, static double complex interpolate_phased_bilinear(const double *ref, const double *phases, + const unsigned char *flags, + const char *sym, float hd, float kd, signed int l) { @@ -181,8 +243,8 @@ static double complex interpolate_phased_bilinear(const double *ref, f = kd - (float)k; assert(f >= 0.0); - val1 = interpolate_phased_linear(ref, phases, hd, k, l); - val2 = interpolate_phased_linear(ref, phases, hd, k+1, l); + val1 = interpolate_phased_linear(ref, phases, flags, sym, hd, k, l); + val2 = interpolate_phased_linear(ref, phases, flags, sym, hd, k+1, l); return (1.0-f)*val1 + f*val2; } @@ -190,6 +252,8 @@ static double complex interpolate_phased_bilinear(const double *ref, static double interpolate_phased_intensity(const double *ref, const double *phases, + const unsigned char *flags, + const char *sym, float hd, float kd, float ld) { signed int l; @@ -201,20 +265,22 @@ static double interpolate_phased_intensity(const double *ref, f = ld - (float)l; assert(f >= 0.0); - val1 = interpolate_phased_bilinear(ref, phases, hd, kd, l); - val2 = interpolate_phased_bilinear(ref, phases, hd, kd, l+1); + val1 = interpolate_phased_bilinear(ref, phases, flags, sym, + hd, kd, l); + val2 = interpolate_phased_bilinear(ref, phases, flags, sym, + hd, kd, l+1); return cabs((1.0-f)*val1 + f*val2); } /* Look up the structure factor for the nearest Bragg condition */ -static double molecule_factor(const double *intensities,const double *phases, - struct rvec q, +static double molecule_factor(const double *intensities, const double *phases, + const unsigned char *flags, struct rvec q, double ax, double ay, double az, double bx, double by, double bz, double cx, double cy, double cz, - GradientMethod m) + GradientMethod m, const char *sym) { float hd, kd, ld; signed int h, k, l; @@ -224,6 +290,9 @@ static double molecule_factor(const double *intensities,const double *phases, kd = q.u * bx + q.v * by + q.w * bz; ld = q.u * cx + q.v * cy + q.w * cz; + /* No flags -> flat intensity distribution */ + if ( flags == NULL ) return 1.0e5; + switch ( m ) { case GRADIENT_MOSAIC : h = (signed int)rint(hd); @@ -232,14 +301,14 @@ static double molecule_factor(const double *intensities,const double *phases, if ( abs(h) > INDMAX ) r = 0.0; else if ( abs(k) > INDMAX ) r = 0.0; else if ( abs(l) > INDMAX ) r = 0.0; - else r = lookup_intensity(intensities, h, k, l); + else r = sym_lookup_intensity(intensities, flags, sym, h, k, l); break; case GRADIENT_INTERPOLATE : - r = interpolate_intensity(intensities, hd, kd, ld); + r = interpolate_intensity(intensities, flags, sym, hd, kd, ld); break; case GRADIENT_PHASED : - r = interpolate_phased_intensity(intensities, phases, - hd, kd, ld); + r = interpolate_phased_intensity(intensities, phases, flags, + sym, hd, kd, ld); break; default: ERROR("This gradient method not implemented yet.\n"); @@ -298,7 +367,8 @@ double water_diffraction(struct rvec q, double en, void get_diffraction(struct image *image, int na, int nb, int nc, const double *intensities, const double *phases, - UnitCell *cell, int do_water, GradientMethod m) + const unsigned char *flags, UnitCell *cell, int do_water, + GradientMethod m, const char *sym) { unsigned int xs, ys; double ax, ay, az; @@ -353,10 +423,10 @@ void get_diffraction(struct image *image, int na, int nb, int nc, I_molecule = 1.0e10; } else { I_molecule = molecule_factor(intensities, - phases, q, + phases, flags, q, ax,ay,az, bx,by,bz,cx,cy,cz, - m); + m, sym); } I_lattice = pow(f_lattice, 2.0); diff --git a/src/diffraction.h b/src/diffraction.h index 5a88bff7..607b066e 100644 --- a/src/diffraction.h +++ b/src/diffraction.h @@ -26,7 +26,8 @@ typedef enum { } GradientMethod; extern void get_diffraction(struct image *image, int na, int nb, int nc, - const double *intensities,const double *phases, - UnitCell *cell, int do_water, GradientMethod m); + const double *intensities, const double *phases, + const unsigned char *flags, UnitCell *cell, + int do_water, GradientMethod m, const char *sym); #endif /* DIFFRACTION_H */ diff --git a/src/indexamajig.c b/src/indexamajig.c index 28573f59..e99b1bd3 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -70,6 +70,8 @@ struct static_index_args IndexingMethod indm; IndexingPrivate *ipriv; const double *intensities; + const unsigned char *flags; + const char *sym; /* Symmetry of "intensities" and "flags" */ struct gpu_context *gctx; int gpu_dev; int peaks; @@ -191,6 +193,7 @@ static void show_help(const char *s) "\nIf you used --simulate, you may also want:\n\n" " --intensities=<file> Specify file containing reflection intensities\n" " to use when simulating.\n" +" -y, --symmetry=<sym> The symmetry of the intensities file.\n" "\n" "\nOptions for greater performance or verbosity:\n\n" " --verbose Be verbose about indexing.\n" @@ -275,22 +278,23 @@ static struct image *get_simage(struct image *template, int alternate) static void simulate_and_write(struct image *simage, struct gpu_context **gctx, - const double *intensities, UnitCell *cell, - int gpu_dev) + const double *intensities, + const unsigned char *flags, UnitCell *cell, + int gpu_dev, const char *sym) { /* Set up GPU if necessary. * Unfortunately, setup has to go here since until now we don't know * enough about the situation. */ if ( (gctx != NULL) && (*gctx == NULL) ) { - *gctx = setup_gpu(0, simage, intensities, gpu_dev); + *gctx = setup_gpu(0, simage, intensities, flags, gpu_dev); } if ( (gctx != NULL) && (*gctx != NULL) ) { get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell); } else { get_diffraction(simage, 24, 24, 40, - intensities, NULL, cell, 0, - GRADIENT_MOSAIC); + intensities, NULL, flags, cell, 0, + GRADIENT_MOSAIC, sym); } record_image(simage, 0); @@ -321,7 +325,9 @@ static void process_image(void *pp, int cookie) int config_polar = pargs->static_args.config_polar; IndexingMethod indm = pargs->static_args.indm; const double *intensities = pargs->static_args.intensities; + const unsigned char *flags = pargs->static_args.flags; struct gpu_context *gctx = pargs->static_args.gctx; + const char *sym = pargs->static_args.sym; image.features = NULL; image.data = NULL; @@ -428,13 +434,13 @@ static void process_image(void *pp, int cookie) if ( config_simulate ) { if ( config_gpu ) { pthread_mutex_lock(pargs->static_args.gpu_mutex); - simulate_and_write(simage, &gctx, intensities, + simulate_and_write(simage, &gctx, intensities, flags, image.indexed_cell, - pargs->static_args.gpu_dev); + pargs->static_args.gpu_dev, sym); pthread_mutex_unlock(pargs->static_args.gpu_mutex); } else { - simulate_and_write(simage, NULL, intensities, - image.indexed_cell, 0); + simulate_and_write(simage, NULL, intensities, flags, + image.indexed_cell, 0, sym); } } @@ -539,6 +545,7 @@ int main(int argc, char *argv[]) char *indm_str = NULL; UnitCell *cell; double *intensities = NULL; + unsigned char *flags; char *intfile = NULL; char *pdb = NULL; char *prefix = NULL; @@ -557,6 +564,7 @@ int main(int argc, char *argv[]) struct beam_params *beam = NULL; double nominal_photon_energy; int gpu_dev = -1; + char *sym = NULL; /* Long options */ const struct option longopts[] = { @@ -579,6 +587,7 @@ int main(int argc, char *argv[]) {"verbose", 0, &config_verbose, 1}, {"alternate", 0, &config_alternate, 1}, {"intensities", 1, NULL, 'q'}, + {"symmetry", 1, NULL, 'y'}, {"pdb", 1, NULL, 'p'}, {"prefix", 1, NULL, 'x'}, {"unpolarized", 0, &config_polar, 0}, @@ -595,7 +604,7 @@ int main(int argc, char *argv[]) }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:o:b:", + while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:o:b:y:", longopts, NULL)) != -1) { switch (c) { @@ -639,6 +648,10 @@ int main(int argc, char *argv[]) threshold = strtof(optarg, NULL); break; + case 'y' : + sym = strdup(optarg); + break; + case 'b' : beam = get_beam_parameters(optarg); if ( beam == NULL ) { @@ -701,6 +714,8 @@ int main(int argc, char *argv[]) } free(outfile); + if ( sym == NULL ) sym = strdup("1"); + if ( speaks == NULL ) { speaks = strdup("zaef"); STATUS("You didn't specify a peak detection method.\n"); @@ -717,12 +732,26 @@ int main(int argc, char *argv[]) free(speaks); if ( intfile != NULL ) { + ReflItemList *items; + int i; + items = read_reflections(intfile, intensities, NULL, NULL, NULL); + + flags = new_list_flag(); + for ( i=0; i<num_items(items); i++ ) { + struct refl_item *it = get_item(items, i); + set_flag(flags, it->h, it->k, it->l, 1); + } + delete_items(items); + } else { + intensities = NULL; + flags = NULL; + } if ( pdb == NULL ) { @@ -859,6 +888,8 @@ int main(int argc, char *argv[]) qargs.static_args.indm = indm; qargs.static_args.ipriv = ipriv; qargs.static_args.intensities = intensities; + qargs.static_args.flags = flags; + qargs.static_args.sym = sym; qargs.static_args.gctx = gctx; qargs.static_args.gpu_dev = gpu_dev; qargs.static_args.peaks = peaks; @@ -881,6 +912,7 @@ int main(int argc, char *argv[]) free(det); cell_free(cell); if ( fh != stdout ) fclose(fh); + free(sym); STATUS("There were %i images. %i could be indexed, of which %i" " looked sane.\n", n_images, qargs.n_indexable, qargs.n_sane); diff --git a/src/pattern_sim.c b/src/pattern_sim.c index b6fcd4e9..41dd3a26 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -68,6 +68,7 @@ static void show_help(const char *s) " this invocation as the given filename.\n" " -i, --intensities=<file> Specify file containing reflection intensities\n" " (and phases) to use.\n" +" -y, --symmetry=<sym> The symmetry of the intensities file.\n" " -t, --gradients=<method> Use <method> for the calculation of shape\n" " transform intensities. Choose from:\n" " mosaic : Take the intensity of the nearest\n" @@ -201,6 +202,7 @@ int main(int argc, char *argv[]) double *intensities; char *rval; double *phases; + unsigned char *flags; int config_simdetails = 0; int config_nearbragg = 0; int config_randomquat = 0; @@ -227,6 +229,7 @@ int main(int argc, char *argv[]) int random_size = 0; double min_size = 0.0; double max_size = 0.0; + char *sym = NULL; /* Long options */ const struct option longopts[] = { @@ -240,6 +243,7 @@ int main(int argc, char *argv[]) {"no-water", 0, &config_nowater, 1}, {"no-noise", 0, &config_nonoise, 1}, {"intensities", 1, NULL, 'i'}, + {"symmetry", 1, NULL, 'y'}, {"powder", 1, NULL, 'w'}, {"gradients", 1, NULL, 't'}, {"pdb", 1, NULL, 'p'}, @@ -254,7 +258,7 @@ int main(int argc, char *argv[]) }; /* Short options */ - while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:g:b:", + while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:g:b:y:", longopts, NULL)) != -1) { switch (c) { @@ -302,6 +306,10 @@ int main(int argc, char *argv[]) beamfile = strdup(optarg); break; + case 'y' : + sym = strdup(optarg); + break; + case 2 : gpu_dev = atoi(optarg); break; @@ -359,6 +367,8 @@ int main(int argc, char *argv[]) } } + if ( sym == NULL ) sym = strdup("1"); + if ( config_simdetails ) { show_details(); return 0; @@ -404,14 +414,20 @@ int main(int argc, char *argv[]) } if ( intfile == NULL ) { + /* Gentle reminder */ STATUS("You didn't specify the file containing the "); STATUS("reflection intensities (with --intensities).\n"); STATUS("I'll simulate a flat intensity distribution.\n"); intensities = NULL; phases = NULL; + flags = NULL; + } else { + + int i; ReflItemList *items; + if ( grad == GRADIENT_PHASED ) { phases = new_list_phase(); } else { @@ -419,9 +435,16 @@ int main(int argc, char *argv[]) } intensities = new_list_intensity(); phases = new_list_phase(); + flags = new_list_flag(); items = read_reflections(intfile, intensities, phases, NULL, NULL); free(intfile); + + for ( i=0; i<num_items(items); i++ ) { + struct refl_item *it = get_item(items, i); + set_flag(flags, it->h, it->k, it->l, 1); + } + delete_items(items); } @@ -516,12 +539,13 @@ int main(int argc, char *argv[]) if ( config_gpu ) { if ( gctx == NULL ) { gctx = setup_gpu(config_nosfac, &image, - intensities, gpu_dev); + intensities, flags, gpu_dev); } get_diffraction_gpu(gctx, &image, na, nb, nc, cell); } else { get_diffraction(&image, na, nb, nc, intensities, phases, - cell, !config_nowater, grad); + flags, cell, !config_nowater, grad, + sym); } if ( image.data == NULL ) { ERROR("Diffraction calculation failed.\n"); @@ -599,6 +623,7 @@ skip: free(intensities); free(outfile); free(filename); + free(sym); return 0; } diff --git a/src/utils.h b/src/utils.h index a8461d6a..1fbd5e4e 100644 --- a/src/utils.h +++ b/src/utils.h @@ -180,6 +180,11 @@ static inline double angle_between(double x1, double y1, double z1, #define TYPE double #include "list_tmp.h" +/* As above, but for simple flags */ +#define LABEL(x) x##_flag +#define TYPE unsigned char +#include "list_tmp.h" + /* ----------- Reflection lists indexed by sequence (not indices) ----------- */ |