Use symmetry when simulating (on the CPU only)

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) ----------- */