Don't try to render PDBs, part I

7 files changed, 109 insertions, 103 deletions

diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 9b2f51f5..cc7b9e5b 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -38,7 +38,7 @@ struct gpu_context
 	cl_command_queue cq;
 	cl_program prog;
 	cl_kernel kern;
-	cl_mem sfacs;
+	cl_mem intensities;
 
 	cl_mem tt;
 	size_t tt_size;
@@ -183,7 +183,7 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 		ERROR("Couldn't set arg 8: %s\n", clError(err));
 		return;
 	}
-	clSetKernelArg(gctx->kern, 9, sizeof(cl_mem), &gctx->sfacs);
+	clSetKernelArg(gctx->kern, 9, sizeof(cl_mem), &gctx->intensities);
 	if ( err != CL_SUCCESS ) {
 		ERROR("Couldn't set arg 9: %s\n", clError(err));
 		return;
@@ -344,7 +344,7 @@ 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,
-                              struct molecule *molecule)
+                              double *intensities)
 {
 	struct gpu_context *gctx;
 	cl_uint nplat;
@@ -352,21 +352,11 @@ struct gpu_context *setup_gpu(int no_sfac, struct image *image,
 	cl_context_properties prop[3];
 	cl_int err;
 	cl_device_id dev;
-	size_t sfac_size;
-	float *sfac_ptr;
+	size_t intensities_size;
+	float *intensities_ptr;
 	size_t maxwgsize;
 	int i;
 
-	if ( molecule == NULL ) return NULL;
-
-	/* Generate structure factors if required */
-	if ( !no_sfac ) {
-		if ( molecule->reflections == NULL ) {
-			get_reflections_cached(molecule,
-			                       ph_lambda_to_en(image->lambda));
-		}
-	}
-
 	STATUS("Setting up GPU..."); fflush(stderr);
 
 	err = clGetPlatformIDs(8, platforms, &nplat);
@@ -411,28 +401,26 @@ struct gpu_context *setup_gpu(int no_sfac, struct image *image,
 	}
 
 	/* Create a single-precision version of the scattering factors */
-	sfac_size = IDIM*IDIM*IDIM*sizeof(cl_float)*2; /* complex */
-	sfac_ptr = malloc(sfac_size);
-	if ( !no_sfac ) {
+	intensities_size = IDIM*IDIM*IDIM*sizeof(cl_float);
+	intensities_ptr = malloc(intensities_size);
+	if ( intensities != NULL ) {
 		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-			sfac_ptr[2*i+0] = creal(molecule->reflections[i]);
-			sfac_ptr[2*i+1] = cimag(molecule->reflections[i]);
+			intensities_ptr[i] = intensities[i];
 		}
 	} else {
 		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-			sfac_ptr[2*i+0] = 10000.0;
-			sfac_ptr[2*i+1] = 0.0;
+			intensities_ptr[i] = 10000.0;
 		}
 	}
-	gctx->sfacs = clCreateBuffer(gctx->ctx,
+	gctx->intensities = clCreateBuffer(gctx->ctx,
 	                             CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
-	                             sfac_size, sfac_ptr, &err);
+	                             intensities_size, intensities_ptr, &err);
 	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't allocate sfac memory\n");
+		ERROR("Couldn't allocate intensities memory\n");
 		free(gctx);
 		return NULL;
 	}
-	free(sfac_ptr);
+	free(intensities_ptr);
 
 	gctx->tt_size = image->width*image->height*sizeof(cl_float);
 	gctx->tt = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, gctx->tt_size,
@@ -478,7 +466,7 @@ void cleanup_gpu(struct gpu_context *gctx)
 	clReleaseProgram(gctx->prog);
 	clReleaseMemObject(gctx->diff);
 	clReleaseMemObject(gctx->tt);
-	clReleaseMemObject(gctx->sfacs);
+	clReleaseMemObject(gctx->intensities);
 
 	/* Release LUTs */
 	for ( i=1; i<=gctx->max_sinc_lut; i++ ) {
diff --git a/src/diffraction-gpu.h b/src/diffraction-gpu.h
index f3bed8bf..b23dab63 100644
--- a/src/diffraction-gpu.h
+++ b/src/diffraction-gpu.h
@@ -26,7 +26,7 @@ struct gpu_context;
 extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
                                 int na, int nb, int nc);
 extern struct gpu_context *setup_gpu(int no_sfac, struct image *image,
-                                     struct molecule *molecule);
+                                     double *intensities);
 extern void cleanup_gpu(struct gpu_context *gctx);
 
 #else
@@ -39,7 +39,7 @@ static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 }
 
 static struct gpu_context *setup_gpu(int no_sfac, struct image *image,
-                                     struct molecule *molecule)
+                                     double *intensities)
 {
 	return NULL;
 }
diff --git a/src/diffraction.c b/src/diffraction.c
index 5938c307..01624be8 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -70,14 +70,14 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az,
 
 
 /* Look up the structure factor for the nearest Bragg condition */
-static double complex molecule_factor(struct molecule *mol, struct rvec q,
+static double molecule_factor(double *intensities, struct rvec q,
                                       double ax, double ay, double az,
                                       double bx, double by, double bz,
                                       double cx, double cy, double cz)
 {
 	double hd, kd, ld;
 	signed int h, k, l;
-	double complex r;
+	double r;
 
 	hd = q.u * ax + q.v * ay + q.w * az;
 	kd = q.u * bx + q.v * by + q.w * bz;
@@ -86,7 +86,7 @@ static double complex molecule_factor(struct molecule *mol, struct rvec q,
 	k = (signed int)rint(kd);
 	l = (signed int)rint(ld);
 
-	r = lookup_sfac(mol->reflections, h, k, l);
+	r = lookup_intensity(intensities, h, k, l);
 
 	return r;
 }
@@ -170,8 +170,8 @@ struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
 }
 
 
-void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac,
-                     int do_water)
+void get_diffraction(struct image *image, int na, int nb, int nc,
+                     double *intensities, int do_water)
 {
 	unsigned int xs, ys;
 	double ax, ay, az;
@@ -188,13 +188,6 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac,
 	/* Allocate (and zero) the "diffraction array" */
 	image->data = calloc(image->width * image->height, sizeof(float));
 
-	if ( !no_sfac ) {
-		if ( image->molecule->reflections == NULL ) {
-			get_reflections_cached(image->molecule,
-			                       ph_lambda_to_en(image->lambda));
-		}
-	}
-
 	/* Needed later for Lorentz calculation */
 	image->twotheta = malloc(image->width * image->height * sizeof(double));
 
@@ -205,8 +198,6 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac,
 	for ( xs=0; xs<image->width*SAMPLING; xs++ ) {
 	for ( ys=0; ys<image->height*SAMPLING; ys++ ) {
 
-		double f_lattice;
-		double complex f_molecule;
 		struct rvec q;
 		float twotheta;
 		double sw = 1.0/(SAMPLING*SAMPLING);  /* Sample weight */
@@ -220,8 +211,9 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac,
 
 			float k;
 			double kw = 1.0/BWSAMPLING;
-			double complex val;
 			double intensity;
+			double f_lattice, I_lattice;
+			double I_molecule;
 
 			/* Calculate k this time round */
 			k = klow + kstep * bwstep;
@@ -233,15 +225,17 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac,
 			                              bx, by, bz,
 			                              cx, cy, cz,
 			                              na, nb, nc);
-			if ( no_sfac ) {
-				f_molecule = 10000.0;
+
+			if ( intensities == NULL ) {
+				I_molecule = 10000.0;
 			} else {
-				f_molecule = molecule_factor(image->molecule, q,
+				I_molecule = molecule_factor(intensities, q,
 			                            ax,ay,az,bx,by,bz,cx,cy,cz);
 			}
 
-			val = f_molecule * f_lattice;
-			intensity = sw * kw * pow(cabs(val), 2.0);
+			I_lattice = pow(f_lattice, 2.0);
+
+			intensity = sw * kw * I_lattice * I_molecule;
 			image->data[x + image->width*y] += intensity;
 
 		}
diff --git a/src/diffraction.h b/src/diffraction.h
index 503e9d19..dd0795fb 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -20,7 +20,8 @@
 #include "cell.h"
 
 extern void get_diffraction(struct image *image, int na, int nb, int nc,
-                            int nosfac, int do_water);
+                            double *intensities, int do_water);
 extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
                          unsigned int sampling, float *ttp, float k);
+
 #endif	/* DIFFRACTION_H */
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 0240826c..29ff0786 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -110,6 +110,8 @@ int main(int argc, char *argv[])
 	char *output = NULL;
 	unsigned int *counts;
 	signed int h, k, l;
+	FILE *fh;
+	char *rval;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -154,7 +156,7 @@ int main(int argc, char *argv[])
 
 	mol = load_molecule();
 	get_reflections_cached(mol, eV_to_J(1790.0));
-	ideal_ref = ideal_intensities(mol->reflections);
+
 	counts = new_list_count();
 
 	if ( template != NULL ) {
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 05cf57d4..4fd7da66 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -30,6 +30,7 @@
 #include "detector.h"
 #include "intensities.h"
 #include "sfac.h"
+#include "reflections.h"
 
 
 static void show_help(const char *s)
@@ -50,6 +51,8 @@ static void show_help(const char *s)
 "                            (a new orientation will be used for each image).\n"
 "     --powder              Write a summed pattern of all images simulated by\n"
 "                            this invocation to results/integr.h5.\n"
+" -i, --intensities=<file>  Specify file containing reflection intensities\n"
+"                            to use.\n"
 "\n"
 "By default, the simulation aims to be as accurate as possible.  For greater\n"
 "speed, or for testing, you can choose to disable certain things using the\n"
@@ -57,7 +60,6 @@ static void show_help(const char *s)
 "\n"
 "     --no-water            Do not simulate water background.\n"
 "     --no-noise            Do not calculate Poisson noise.\n"
-"     --no-sfac             Pretend that all structure factors are 1.\n"
 );
 }
 
@@ -78,21 +80,15 @@ static void show_details()
 "na = number of unit cells in 'a' direction (likewise nb, nc)\n"
 " q = reciprocal vector (1/d convention, not 2pi/d)\n"
 "\n"
-"This value is multiplied by the complex structure factor at the nearest\n"
-"Bragg position, i.e. the gradient of the shape transform across each\n"
-"appearance of the shape transform is not included, for speed of calculation.\n"
+"This square modulus of this value is take, and multiplied by the Bragg\n"
+"intensitiy at the neares Bragg position.  That means that the gradient of\n"
+"the underlying molecule transform is not included, for speed of calculation.\n"
+"The Bragg intensities are taken from the file specified on the command line\n"
+"with the --intensities option.\n"
 "\n"
-"Complex structure factors are calculated using a combination of the Henke\n"
-"and Waasmeier-Kirfel scattering factors. The Henke factors are complex\n"
-"and energy dependence, whereas the Waas-Kirf values are real-valued and\n"
-"|q|-dependent.  The difference between the Waas-Kirf value at the\n"
-"appropriate |q| and the same value at |q|=0 is subtracted from the Henke\n"
-"value.  The Henke values are linearly interpolated from the provided tables\n"
-"(note that the interpolation should really be exponential).\n"
-"\n"
-"The modulus of the structure factor is taken and squared.  Intensity from\n"
-"water is then added according to the first term of equation 5 from\n"
-"Phys Chem Chem Phys 2003 (5) 1981--1991.\n"
+"Intensity from water is then added according to the first term of equation\n"
+"5 from Phys Chem Chem Phys 2003 (5) 1981--1991.  This simulates the\n"
+"coherent, elastic part of the diffuse scattering from the water jet only.\n"
 "\n"
 "Expected intensities at the CCD are then calculated using:\n"
 "\n"
@@ -101,14 +97,12 @@ static void show_details()
 "I0 = number of photons per unit area in the incident beam\n"
 " r = Thomson radius\n"
 " S = solid angle of corresponding pixel\n"
+"where |F(q)|^2 is the value calculated as described above.\n"
 "\n"
 "Poisson counts are generated from the expected intensities using Knuth's\n"
 "algorithm.  When the intensity is sufficiently high that Knuth's algorithm\n"
 "would result in machine precision problems, a normal distribution with\n"
 "standard deviation sqrt(I) is used instead.\n"
-"\n"
-"The coherent, elastic part of the diffuse scattering from the water jet can\n"
-"be simulated.\n"
 );
 }
 
@@ -153,6 +147,8 @@ int main(int argc, char *argv[])
 	struct image image;
 	struct gpu_context *gctx = NULL;
 	double *powder;
+	char *intfile = NULL;
+	double *intensities;
 	int config_simdetails = 0;
 	int config_nearbragg = 0;
 	int config_randomquat = 0;
@@ -178,7 +174,7 @@ int main(int argc, char *argv[])
 		{"no-images",          0, &config_noimages,    1},
 		{"no-water",           0, &config_nowater,     1},
 		{"no-noise",           0, &config_nonoise,     1},
-		{"no-sfac",            0, &config_nosfac,      1},
+		{"intensities",        0, NULL,               'i'},
 		{"powder",             0, &config_powder,      1},
 		{0, 0, NULL, 0}
 	};
@@ -202,6 +198,11 @@ int main(int argc, char *argv[])
 			break;
 		}
 
+		case 'i' : {
+			intfile = strdup(optarg);
+			break;
+		}
+
 		case 0 : {
 			break;
 		}
@@ -224,6 +225,17 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	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;
+	} else {
+		intensities = read_reflections(intfile);
+		free(intfile);
+	}
+
 	/* Define image parameters */
 	image.width = 1024;
 	image.height = 1024;
@@ -276,11 +288,11 @@ int main(int argc, char *argv[])
 		if ( config_gpu ) {
 			if ( gctx == NULL ) {
 				gctx = setup_gpu(config_nosfac, &image,
-				                 image.molecule);
+				                 intensities);
 			}
 			get_diffraction_gpu(gctx, &image, na, nb, nc);
 		} else {
-			get_diffraction(&image, na, nb, nc, config_nosfac,
+			get_diffraction(&image, na, nb, nc, intensities,
 			                !config_nowater);
 		}
 		if ( image.molecule == NULL ) {
diff --git a/src/process_hkl.c b/src/process_hkl.c
index d110695b..53ac1169 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -37,31 +37,31 @@ static void show_help(const char *s)
 	printf(
 "Assemble and process FEL Bragg intensities.\n"
 "\n"
-"  -h, --help              Display this help message.\n"
-"  -i, --input=<filename>  Specify input filename (\"-\" for stdin).\n"
-"  -o, --output=<filename> Specify output filename for merged intensities\n"
-"                           (don't specify for no output).\n"
+"  -h, --help                Display this help message.\n"
+"  -i, --input=<filename>    Specify input filename (\"-\" for stdin).\n"
+"  -o, --output=<filename>   Specify output filename for merged intensities\n"
+"                             (don't specify for no output).\n"
 "\n"
-"      --max-only          Take the integrated intensity to be equal to the\n"
-"                           maximum intensity measured for that reflection.\n"
-"                           The default is to use the mean value from all\n"
-"                           measurements.\n"
-"  -e, --output-every=<n>  Analyse figures of merit after every n patterns\n"
-"                           Default: 1000.  A value of zero means to do the\n"
-"                           analysis only after reading all the patterns.\n"
-"      --no-analyse        Don't perform any kind of analysis, just merge the\n"
-"                           intensities.\n"
-"      --sum               Sum (rather than average) the intensities for the\n"
-"                           final output list.  This is useful for comparing\n"
-"                           results to radially summed powder patterns, but\n"
-"                           will break R-factor analysis.\n"
-"  -r, --rvsq              Output lists of R vs |q| (\"Luzzatti plots\") when\n"
-"                           analysing figures of merit.\n"
-"      --stop-after=<n>    Stop after processing n patterns.  Zero means\n"
-"                           keep going until the end of the input, and is the\n"
-"                           default.\n"
-"      --zone-axis         Output an [001] zone axis pattern each time the\n"
-"                           figures of merit are analysed.\n");
+"      --max-only            Take the integrated intensity to be equal to the\n"
+"                             maximum intensity measured for that reflection.\n"
+"                             The default is to use the mean value from all\n"
+"                             measurements.\n"
+"      --sum                 Sum (rather than average) the intensities for the\n"
+"                             final output list.  This is useful for comparing\n"
+"                             results to radially summed powder patterns, but\n"
+"                             will break R-factor analysis.\n"
+"      --stop-after=<n>      Stop after processing n patterns.  Zero means\n"
+"                             keep going until the end of the input, and is\n"
+"                             the default.\n"
+"  -c, --compare-with=<file> Compare with reflection intensities in this file\n"
+"\n"
+"  -e, --output-every=<n>    Analyse figures of merit after every n patterns\n"
+"                             Default: 1000.  A value of zero means to do the\n"
+"                             analysis only after reading all the patterns.\n"
+"  -r, --rvsq                Output lists of R vs |q| (\"Luzzatti plots\")\n"
+"                             when analysing figures of merit.\n"
+"      --zone-axis           Output an [001] zone axis pattern each time the\n"
+"                             figures of merit are analysed.\n");
 }
 
 
@@ -193,8 +193,8 @@ int main(int argc, char *argv[])
 	int config_rvsq = 0;
 	int config_stopafter = 0;
 	int config_zoneaxis = 0;
-	int config_noanalyse = 0;
 	int config_sum = 0;
+	char *intfile = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -203,10 +203,10 @@ int main(int argc, char *argv[])
 		{"output",             1, NULL,               'o'},
 		{"max-only",           0, &config_maxonly,     1},
 		{"output-every",       1, NULL,               'e'},
-		{"no-analyse",         0, &config_noanalyse,   1},
 		{"rvsq",               0, NULL,               'r'},
 		{"stop-after",         1, NULL,               's'},
 		{"zone-axis",          0, &config_zoneaxis,    1},
+		{"compare-with",       0, NULL,               'c'},
 		{"sum",                0, &config_sum,    1},
 		{0, 0, NULL, 0}
 	};
@@ -245,6 +245,11 @@ int main(int argc, char *argv[])
 			break;
 		}
 
+		case 'c' : {
+			intfile = strdup(optarg);
+			break;
+		}
+
 		case 0 : {
 			break;
 		}
@@ -261,14 +266,18 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	if ( intfile != NULL ) {
+		STATUS("Comparing against '%s'\n", intfile);
+		trueref = read_reflections(intfile);
+		free(intfile);
+	} else {
+		trueref = NULL;
+	}
+
 	ref = new_list_intensity();
 	counts = new_list_count();
 
 	mol = load_molecule();
-	if ( !config_noanalyse ) {
-		get_reflections_cached(mol, eV_to_J(2.0e3));
-		trueref = ideal_intensities(mol->reflections);
-	}
 
 	if ( strcmp(filename, "-") == 0 ) {
 		fh = stdin;
@@ -331,7 +340,7 @@ int main(int argc, char *argv[])
 
 	fclose(fh);
 
-	if ( !config_noanalyse ) {
+	if ( trueref != NULL ) {
 		process_reflections(ref, trueref, counts, n_patterns, mol->cell,
 		                    config_rvsq, config_zoneaxis);
 	}