pattern_sim.c 11.9 KB
Newer Older
1
/*
Thomas White's avatar
Thomas White committed
2
3
4
 * pattern_sim.c
 *
 * Simulate diffraction patterns from small crystals
5
 *
Thomas White's avatar
Thomas White committed
6
 * (c) 2006-2010 Thomas White <taw@physics.org>
7
 *
Thomas White's avatar
Thomas White committed
8
 * Part of CrystFEL - crystallography with a FEL
9
10
11
12
13
14
15
16
17
18
19
20
21
 *
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
Thomas White's avatar
Thomas White committed
22
#include <getopt.h>
23

24
#include "image.h"
Thomas White's avatar
Thomas White committed
25
#include "diffraction.h"
26
#include "diffraction-gpu.h"
27
#include "cell.h"
28
29
#include "utils.h"
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
30
#include "detector.h"
31
#include "peaks.h"
Thomas White's avatar
Thomas White committed
32
#include "sfac.h"
33
#include "reflections.h"
34
#include "parameters-lcls.tmp"
35
36


Thomas White's avatar
Thomas White committed
37
static void show_help(const char *s)
38
{
39
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
40
	printf(
Thomas White's avatar
Thomas White committed
41
"Simulate diffraction patterns from small crystals probed with femtosecond\n"
Thomas White's avatar
Thomas White committed
42
43
"pulses of X-rays from a free electron laser.\n"
"\n"
44
" -h, --help                Display this help message.\n"
45
46
47
" -p, --pdb=<file>          PDB file from which to get the unit cell.\n"
"                            (The actual Bragg intensities come from the\n"
"                            intensities file)\n"
48
"     --simulation-details  Show technical details of the simulation.\n"
49
"     --gpu                 Use the GPU to speed up the calculation.\n"
50
51
52
53
"\n"
"     --near-bragg          Output h,k,l,I near Bragg conditions.\n"
" -n, --number=<N>          Generate N images.  Default 1.\n"
"     --no-images           Do not output any HDF5 files.\n"
54
55
56
57
58
59
" -o, --output=<filename>   Output HDF5 filename.  Default: sim.h5.\n"
"                            If you choose to simulate more than one pattern,\n"
"                            the filename given will be postfixed with a\n"
"                            hyphen, the image number and '.h5'.  In this\n"
"                            case, the default value is 'sim', such that the\n"
"                            files produced are sim-1.h5, sim-2.h5 and so on.\n"
Thomas White's avatar
Thomas White committed
60
" -r, --random-orientation  Use a randomly generated orientation\n"
61
"                            (a new orientation will be used for each image).\n"
Thomas White's avatar
Thomas White committed
62
63
"     --powder              Write a summed pattern of all images simulated by\n"
"                            this invocation to results/integr.h5.\n"
64
65
" -i, --intensities=<file>  Specify file containing reflection intensities\n"
"                            to use.\n"
66
67
68
69
70
71
72
73
" -g, --gradients=<method>  Use <method> for the calculation of shape\n"
"                            transform intensities.  Choose from:\n"
"                             mosaic      : Take the intensity of the nearest\n"
"                                           Bragg position.\n"
"                             interpolate : Interpolate trilinearly between\n"
"                                           six adjacent Bragg intensities.\n"
"                             phased      : As 'interpolate', but take phase\n"
"                                           values into account.\n"
Thomas White's avatar
Thomas White committed
74
75
76
77
78
"\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"
"following options.\n"
"\n"
79
"     --no-water            Do not simulate water background.\n"
Thomas White's avatar
Thomas White committed
80
"     --no-noise            Do not calculate Poisson noise.\n"
81
);
Thomas White's avatar
Thomas White committed
82
83
84
}


85
static void show_details()
Thomas White's avatar
Thomas White committed
86
{
Thomas White's avatar
Thomas White committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
	printf(
"This program simulates diffraction patterns from small crystals illuminated\n"
"with femtosecond X-ray pulses from a free electron laser.\n"
"\n"
"The lattice transform from the specified number of unit cells is calculated\n"
"using the closed-form solution for a truncated lattice:\n"
"\n"
"F(q) =  sin(pi*na*q.a)/sin(pi*q.a)\n"
"      * sin(pi*nb*q.b)/sin(pi*q.b)\n"
"      * sin(pi*nc*q.c)/sin(pi*q.c)\n"
"\n"
"na = number of unit cells in 'a' direction (likewise nb, nc)\n"
" q = reciprocal vector (1/d convention, not 2pi/d)\n"
"\n"
101
"This square modulus of this value is taken, and multiplied by the Bragg\n"
102
103
104
105
"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"
Thomas White's avatar
Thomas White committed
106
"\n"
107
108
109
"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"
Thomas White's avatar
Thomas White committed
110
111
112
113
114
115
116
117
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
"I(q) = I0 * r^2 * |F(q)|^2 * S\n"
"\n"
"I0 = number of photons per unit area in the incident beam\n"
" r = Thomson radius\n"
" S = solid angle of corresponding pixel\n"
118
"where |F(q)|^2 is the value calculated as described above.\n"
Thomas White's avatar
Thomas White committed
119
120
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
121
122
123
"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"
124
);
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
}


static struct quaternion read_quaternion()
{
	do {

		int r;
		float w, x, y, z;
		char line[1024];
		char *rval;

		printf("Please input quaternion: w x y z\n");
		rval = fgets(line, 1023, stdin);
		if ( rval == NULL ) return invalid_quaternion();
		chomp(line);

		r = sscanf(line, "%f %f %f %f", &w, &x, &y, &z);
		if ( r == 4 ) {

			struct quaternion quat;

			quat.w = w;
			quat.x = x;
			quat.y = y;
			quat.z = z;

			return quat;

		} else {
155
			ERROR("Invalid rotation '%s'\n", line);
156
157
158
		}

	} while ( 1 );
159
160
161
162
163
}


int main(int argc, char *argv[])
{
164
	int c;
165
	struct image image;
166
	struct gpu_context *gctx = NULL;
167
	double *powder;
168
169
	char *intfile = NULL;
	double *intensities;
170
	double *phases;
Thomas White's avatar
Thomas White committed
171
	int config_simdetails = 0;
172
173
	int config_nearbragg = 0;
	int config_randomquat = 0;
174
	int config_noimages = 0;
175
	int config_nowater = 0;
Thomas White's avatar
Thomas White committed
176
	int config_nonoise = 0;
177
	int config_nosfac = 0;
178
	int config_gpu = 0;
179
	int config_powder = 0;
180
	char *filename = NULL;
181
	char *grad_str = NULL;
182
	char *outfile = NULL;
183
	GradientMethod grad;
Thomas White's avatar
Thomas White committed
184
185
186
	int ndone = 0;    /* Number of simulations done (images or not) */
	int number = 1;   /* Number used for filename of image */
	int n_images = 1; /* Generate one image by default */
187
	int done = 0;
188
	UnitCell *cell;
189
	unsigned int *counts;
Thomas White's avatar
Thomas White committed
190

191
	/* Long options */
Thomas White's avatar
Thomas White committed
192
	const struct option longopts[] = {
193
194
		{"help",               0, NULL,               'h'},
		{"simulation-details", 0, &config_simdetails,  1},
195
		{"gpu",                0, &config_gpu,         1},
196
197
198
		{"near-bragg",         0, &config_nearbragg,   1},
		{"random-orientation", 0, NULL,               'r'},
		{"number",             1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
199
		{"no-images",          0, &config_noimages,    1},
200
		{"no-water",           0, &config_nowater,     1},
Thomas White's avatar
Thomas White committed
201
		{"no-noise",           0, &config_nonoise,     1},
202
		{"intensities",        1, NULL,               'i'},
203
		{"powder",             0, &config_powder,      1},
204
		{"gradients",          1, NULL,               'g'},
205
		{"pdb",                1, NULL,               'p'},
206
		{"output",             1, NULL,               'o'},
207
		{0, 0, NULL, 0}
Thomas White's avatar
Thomas White committed
208
	};
209

210
	/* Short options */
211
	while ((c = getopt_long(argc, argv, "hrn:i:g:p:o:", longopts, NULL)) != -1) {
212

Thomas White's avatar
Thomas White committed
213
214
215
216
217
		switch (c) {
		case 'h' : {
			show_help(argv[0]);
			return 0;
		}
218

219
220
221
222
223
224
225
226
227
228
		case 'r' : {
			config_randomquat = 1;
			break;
		}

		case 'n' : {
			n_images = atoi(optarg);
			break;
		}

229
230
231
232
233
		case 'i' : {
			intfile = strdup(optarg);
			break;
		}

234
235
236
237
238
		case 'g' : {
			grad_str = strdup(optarg);
			break;
		}

239
240
241
242
243
		case 'p' : {
			filename = strdup(optarg);
			break;
		}

244
245
246
247
248
		case 'o' : {
			outfile = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
249
250
251
		case 0 : {
			break;
		}
252

Thomas White's avatar
Thomas White committed
253
254
255
		default : {
			return 1;
		}
256
257
258
259
		}

	}

260
261
262
263
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

264
265
266
267
268
269
270
271
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

Thomas White's avatar
Thomas White committed
272
	if ( config_simdetails ) {
273
		show_details();
Thomas White's avatar
Thomas White committed
274
275
276
		return 0;
	}

277
278
279
280
281
282
	if ( (!config_nowater) && config_gpu ) {
		ERROR("Cannot simulate water scattering on the GPU.\n");
		ERROR("Please try again with the --no-water option.\n");
		return 1;
	}

283
284
285
286
287
288
	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;
289
		counts = NULL;
290
		phases = NULL;
291
	} else {
292
		counts = new_list_count();
293
294
		phases = new_list_phase();
		intensities = read_reflections(intfile, counts, phases);
295
296
297
		free(intfile);
	}

298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
	if ( grad_str == NULL ) {
		STATUS("You didn't specify a gradient calculation method, so"
		       " I'm using the 'mosaic' method, which is fastest.\n");
		grad = GRADIENT_MOSAIC;
	} else if ( strcmp(grad_str, "mosaic") == 0 ) {
		grad = GRADIENT_MOSAIC;
	} else if ( strcmp(grad_str, "interpolate") == 0) {
		grad = GRADIENT_INTERPOLATE;
	} else if ( strcmp(grad_str, "phased") == 0) {
		grad = GRADIENT_PHASED;
	} else {
		ERROR("Unrecognised gradient method '%s'\n", grad_str);
		return 1;
	}
	free(grad_str);

	if ( config_gpu && (grad != GRADIENT_MOSAIC) ) {
		ERROR("Only the mosaic method can be used for gradients when"
		      "calculating on the GPU.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
320
	/* Define image parameters */
Thomas White's avatar
Thomas White committed
321
322
	image.width = 1024;
	image.height = 1024;
323
	image.lambda = ph_en_to_lambda(eV_to_J(PHOTON_ENERGY)); /* Wavelength */
324
	cell = load_cell_from_pdb(filename);
Thomas White's avatar
Thomas White committed
325
326
327
	if ( cell == NULL ) {
		exit(1);
	}
328
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
329
330
	image.features = NULL;
	image.flags = NULL;
331
332
	image.f0 = 1.0;
	image.f0_available = 1;
Thomas White's avatar
Thomas White committed
333

Thomas White's avatar
Thomas White committed
334
	#include "geometry-lcls.tmp"
335

336
337
	powder = calloc(image.width*image.height, sizeof(*powder));

Thomas White's avatar
Thomas White committed
338
	/* Splurge a few useful numbers */
339
	STATUS("Wavelength is %f nm\n", image.lambda/1.0e-9);
Thomas White's avatar
Thomas White committed
340

341
342
	do {

Thomas White's avatar
Thomas White committed
343
		int na, nb, nc;
344
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
345

346
347
348
		//na = 8*random()/RAND_MAX + 4;
		//nb = 8*random()/RAND_MAX + 4;
		//nc = 16*random()/RAND_MAX + 30;
Thomas White's avatar
Thomas White committed
349
350
351
		na = 24;
		nb = 24;
		nc = 40;
Thomas White's avatar
Thomas White committed
352

353
354
355
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
			image.orientation = random_quaternion();
356
		} else {
357
			image.orientation = read_quaternion();
358
359
		}

Thomas White's avatar
Thomas White committed
360
361
362
363
		STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
		       image.orientation.w, image.orientation.x,
		       image.orientation.y, image.orientation.z);

Thomas White's avatar
D'oh    
Thomas White committed
364
		if ( !quaternion_valid(image.orientation) ) {
365
			ERROR("Orientation modulus is not zero!\n");
366
367
368
369
370
371
			return 1;
		}

		/* Ensure no residual information */
		image.data = NULL;
		image.twotheta = NULL;
372

373
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
374
375
376
		STATUS("Particle size = %i x %i x %i (=%5.2f x %5.2f x %5.2f nm)\n",
	               na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9);

377
		if ( config_gpu ) {
378
379
			if ( gctx == NULL ) {
				gctx = setup_gpu(config_nosfac, &image,
380
				                 intensities, counts);
381
			}
382
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
383
		} else {
384
			get_diffraction(&image, na, nb, nc, intensities, counts,
385
			                phases, cell, !config_nowater, grad);
386
		}
387
		if ( image.data == NULL ) {
388
389
390
391
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

392
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
393

Thomas White's avatar
Thomas White committed
394
		if ( config_nearbragg ) {
Thomas White's avatar
Thomas White committed
395
			find_projected_peaks(&image, cell);
396
			output_intensities(&image, cell, NULL, 1);
Thomas White's avatar
Thomas White committed
397
398
		}

399
400
401
402
403
404
405
406
		if ( config_powder ) {

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
407
				powder[x+w*y] += (double)image.data[x+w*y];
408
409
410
411
			}
			}

			if ( !(ndone % 10) ) {
412
				hdf5_write("results/integr.h5", powder,
413
				           image.width, image.height,
414
				           H5T_NATIVE_DOUBLE);
415
416
417
			}
		}

418
419
420
421
		if ( !config_noimages ) {

			char filename[1024];

422
423
424
425
426
427
428
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

429
			number++;
Thomas White's avatar
Thomas White committed
430

431
432
			/* Write the output file */
			hdf5_write(filename, image.data,
433
			           image.width, image.height, H5T_NATIVE_FLOAT);
434
435

		}
436

437
438
439
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
440

441
skip:
Thomas White's avatar
Thomas White committed
442
		ndone++;
443

444
		if ( n_images && (ndone >= n_images) ) done = 1;
Thomas White's avatar
Thomas White committed
445

446
	} while ( !done );
447

448
449
450
451
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

452
453
	free(image.det.panels);
	free(powder);
454
455
	free(cell);
	free(intensities);
456

457
	return 0;
458
}