pattern_sim.c 14.3 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 "beam-parameters.h"
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"
Thomas White's avatar
Thomas White committed
45
"\n"
46
47
48
" -p, --pdb=<file>          PDB file from which to get the unit cell.\n"
"                            (The actual Bragg intensities come from the\n"
"                            intensities file)\n"
49
"     --simulation-details  Show technical details of the simulation.\n"
50
"     --gpu                 Use the GPU to speed up the calculation.\n"
51
52
" -g, --geometry=<file>     Get detector geometry from file.\n"
" -b, --beam=<file>         Get beam parameters from file.\n"
53
54
55
56
"\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"
57
58
59
60
61
62
" -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
63
" -r, --random-orientation  Use a randomly generated orientation\n"
64
"                            (a new orientation will be used for each image).\n"
65
66
"     --powder=<file>       Write a summed pattern of all images simulated by\n"
"                            this invocation as the given filename.\n"
67
" -i, --intensities=<file>  Specify file containing reflection intensities\n"
Thomas White's avatar
Thomas White committed
68
"                            (and phases) to use.\n"
Thomas White's avatar
Thomas White committed
69
" -t, --gradients=<method>  Use <method> for the calculation of shape\n"
70
71
"                            transform intensities.  Choose from:\n"
"                             mosaic      : Take the intensity of the nearest\n"
Thomas White's avatar
Thomas White committed
72
73
74
"                                           Bragg position.  This is the\n"
"                                           fastest method and the only one\n"
"                                           supported on the GPU.\n"
75
76
"                             interpolate : Interpolate trilinearly between\n"
"                                           six adjacent Bragg intensities.\n"
Thomas White's avatar
Thomas White committed
77
78
"                                           This method has intermediate\n"
"                                           accuracy.\n"
79
"                             phased      : As 'interpolate', but take phase\n"
Thomas White's avatar
Thomas White committed
80
81
82
"                                           values into account.  This is the\n"
"                                           most accurate method, but the\n"
"                                           slowest.\n"
83
84
85
"     --really-random       Use really random numbers for the orientation and\n"
"                            crystal size.  By default, the same sequence\n"
"                            will be used for each new run.\n"
Thomas White's avatar
Thomas White committed
86
87
88
89
90
"\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"
91
"     --no-water            Do not simulate water background.\n"
Thomas White's avatar
Thomas White committed
92
"     --no-noise            Do not calculate Poisson noise.\n"
93
);
Thomas White's avatar
Thomas White committed
94
95
96
}


97
static void show_details()
Thomas White's avatar
Thomas White committed
98
{
Thomas White's avatar
Thomas White committed
99
100
101
102
103
	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"
104
105
"using the closed-form solution for a truncated lattice faceted on the\n"
"(001), (010) and (100) planes:\n"
Thomas White's avatar
Thomas White committed
106
"\n"
107
108
109
"I_latt(q) =  sin^2(pi*na*q.a)/sin^2(pi*q.a)\n"
"           * sin^2(pi*nb*q.b)/sin^2(pi*q.b)\n"
"           * sin^2(pi*nc*q.c)/sin^2(pi*q.c)\n"
Thomas White's avatar
Thomas White committed
110
111
112
113
"\n"
"na = number of unit cells in 'a' direction (likewise nb, nc)\n"
" q = reciprocal vector (1/d convention, not 2pi/d)\n"
"\n"
114
115
116
"This is multiplied by a model of the underlying molecular transform, I_mol(q).\n"
"This can be approximated to varying levels of accuracy by the methods given by\n"
"the '--gradients' option.\n"
Thomas White's avatar
Thomas White committed
117
"\n"
118
119
"Intensity from water is added according to the first term of equation 5\n"
"from Phys Chem Chem Phys 2003 (5) 1981--1991.  This simulates the\n"
120
"coherent, elastic part of the diffuse scattering from the water jet only.\n"
Thomas White's avatar
Thomas White committed
121
122
123
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
124
"I(q) = I0 * r^2 * I_latt(q) * I_mol(q) * S\n"
Thomas White's avatar
Thomas White committed
125
126
127
128
"\n"
"I0 = number of photons per unit area in the incident beam\n"
" r = Thomson radius\n"
" S = solid angle of corresponding pixel\n"
129
130
131
"\n"
"Polarisation is not currently included in pattern_sim, although it is included\n"
"in the analysis of Bragg peaks by 'indexamajig'.\n"
Thomas White's avatar
Thomas White committed
132
133
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
134
135
136
"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"
137
);
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
}


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 {
168
			ERROR("Invalid rotation '%s'\n", line);
169
170
171
		}

	} while ( 1 );
172
173
174
175
176
}


int main(int argc, char *argv[])
{
177
	int c;
178
	struct image image;
179
	struct gpu_context *gctx = NULL;
180
	double *powder;
181
182
	char *intfile = NULL;
	double *intensities;
183
	double *phases;
Thomas White's avatar
Thomas White committed
184
	int config_simdetails = 0;
185
186
	int config_nearbragg = 0;
	int config_randomquat = 0;
187
	int config_noimages = 0;
188
	int config_nowater = 0;
Thomas White's avatar
Thomas White committed
189
	int config_nonoise = 0;
190
	int config_nosfac = 0;
191
	int config_gpu = 0;
192
	int config_random = 0;
193
	char *powder_fn = NULL;
194
	char *filename = NULL;
195
	char *grad_str = NULL;
196
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
197
	char *geometry = NULL;
198
	char *beamfile = NULL;
199
	GradientMethod grad;
Thomas White's avatar
Thomas White committed
200
201
202
	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 */
203
	int done = 0;
Thomas White's avatar
Thomas White committed
204
205
	UnitCell *input_cell;
	struct quaternion orientation;
Thomas White's avatar
Thomas White committed
206
	int gpu_dev = -1;
Thomas White's avatar
Thomas White committed
207

208
	/* Long options */
Thomas White's avatar
Thomas White committed
209
	const struct option longopts[] = {
210
211
		{"help",               0, NULL,               'h'},
		{"simulation-details", 0, &config_simdetails,  1},
212
		{"gpu",                0, &config_gpu,         1},
213
214
215
		{"near-bragg",         0, &config_nearbragg,   1},
		{"random-orientation", 0, NULL,               'r'},
		{"number",             1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
216
		{"no-images",          0, &config_noimages,    1},
217
		{"no-water",           0, &config_nowater,     1},
Thomas White's avatar
Thomas White committed
218
		{"no-noise",           0, &config_nonoise,     1},
219
		{"intensities",        1, NULL,               'i'},
220
		{"powder",             1, NULL,               'w'},
Thomas White's avatar
Thomas White committed
221
		{"gradients",          1, NULL,               't'},
222
		{"pdb",                1, NULL,               'p'},
223
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
224
		{"geometry",           1, NULL,               'g'},
225
		{"beam",               1, NULL,               'b'},
226
		{"really-random",      0, &config_random,      1},
Thomas White's avatar
Thomas White committed
227
		{"gpu-dev",            1, NULL,                2},
228
		{0, 0, NULL, 0}
Thomas White's avatar
Thomas White committed
229
	};
230

231
	/* Short options */
232
	while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:g:b:",
Thomas White's avatar
Thomas White committed
233
	                        longopts, NULL)) != -1) {
234

Thomas White's avatar
Thomas White committed
235
		switch (c) {
Thomas White's avatar
Thomas White committed
236
		case 'h' :
Thomas White's avatar
Thomas White committed
237
238
			show_help(argv[0]);
			return 0;
239

Thomas White's avatar
Thomas White committed
240
		case 'r' :
241
242
243
			config_randomquat = 1;
			break;

Thomas White's avatar
Thomas White committed
244
		case 'n' :
245
246
247
			n_images = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
248
		case 'i' :
249
250
251
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
252
		case 't' :
253
254
255
			grad_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
256
		case 'p' :
257
258
259
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
260
		case 'o' :
261
262
263
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
264
		case 'w' :
265
266
267
			powder_fn = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
268
		case 'g' :
Thomas White's avatar
Thomas White committed
269
270
271
			geometry = strdup(optarg);
			break;

272
273
274
275
		case 'b' :
			beamfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
276
277
278
279
		case 2 :
			gpu_dev = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
280
		case 0 :
Thomas White's avatar
Thomas White committed
281
			break;
282

Thomas White's avatar
Thomas White committed
283
		default :
Thomas White's avatar
Thomas White committed
284
285
			return 1;
		}
286
287
288

	}

289
	if ( config_random ) {
290
291
292
293
294
295
296
297
		FILE *fh;
		unsigned int seed;
		fh = fopen("/dev/urandom", "r");
		fread(&seed, sizeof(seed), 1, fh);
		fclose(fh);
		srand(seed);
	}

298
299
300
301
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

302
303
304
305
306
307
308
309
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

Thomas White's avatar
Thomas White committed
310
	if ( config_simdetails ) {
311
		show_details();
Thomas White's avatar
Thomas White committed
312
313
314
		return 0;
	}

315
316
317
318
319
320
	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;
	}

321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
	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
343
344
345
346
347
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

348
349
350
351
352
353
	if ( beamfile == NULL ) {
		ERROR("You need to specify a beam parameter file"
		      " with --beam\n");
		return 1;
	}

354
355
356
357
358
359
360
361
	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;
	} else {
Thomas White's avatar
Thomas White committed
362
		ReflItemList *items;
363
364
365
366
367
		if ( grad == GRADIENT_PHASED ) {
			phases = new_list_phase();
		} else {
			phases = NULL;
		}
Thomas White's avatar
Thomas White committed
368
369
		intensities = new_list_intensity();
		phases = new_list_phase();
370
371
		items = read_reflections(intfile, intensities, phases,
		                         NULL, NULL);
372
		free(intfile);
Thomas White's avatar
Thomas White committed
373
		delete_items(items);
374
375
	}

376
377
378
379
380
381
382
	image.det = get_detector_geometry(geometry);
	if ( image.det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

383
384
385
386
387
388
389
	image.beam = get_beam_parameters(beamfile);
	if ( image.beam == NULL ) {
		ERROR("Failed to read beam parameters from '%s'\n", beamfile);
		return 1;
	}
	free(beamfile);

Thomas White's avatar
Thomas White committed
390
	/* Define image parameters */
391
392
	image.width = image.det->max_x + 1;
	image.height = image.det->max_y + 1;
393
	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
Thomas White's avatar
Thomas White committed
394
395
396
397

	/* Load unit cell */
	input_cell = load_cell_from_pdb(filename);
	if ( input_cell == NULL ) {
Thomas White's avatar
Thomas White committed
398
399
		exit(1);
	}
Thomas White's avatar
Thomas White committed
400
401

	/* Initialise stuff */
402
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
403
404
	image.features = NULL;
	image.flags = NULL;
405
406
	image.f0 = 1.0;
	image.f0_available = 1;
Thomas White's avatar
Thomas White committed
407

408
409
	powder = calloc(image.width*image.height, sizeof(*powder));

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

413
414
	do {

Thomas White's avatar
Thomas White committed
415
		int na, nb, nc;
416
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
417
		UnitCell *cell;
Thomas White's avatar
Thomas White committed
418

419
420
421
		//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
422
423
424
		na = 24;
		nb = 24;
		nc = 40;
Thomas White's avatar
Thomas White committed
425

426
427
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
Thomas White's avatar
Thomas White committed
428
			orientation = random_quaternion();
429
		} else {
Thomas White's avatar
Thomas White committed
430
			orientation = read_quaternion();
431
432
		}

Thomas White's avatar
Thomas White committed
433
		STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
Thomas White's avatar
Thomas White committed
434
435
		       orientation.w, orientation.x,
		       orientation.y, orientation.z);
Thomas White's avatar
Thomas White committed
436

Thomas White's avatar
Thomas White committed
437
		if ( !quaternion_valid(orientation) ) {
438
			ERROR("Orientation modulus is not zero!\n");
439
440
441
			return 1;
		}

Thomas White's avatar
Thomas White committed
442
443
		cell = cell_rotate(input_cell, orientation);

444
445
446
		/* Ensure no residual information */
		image.data = NULL;
		image.twotheta = NULL;
447

448
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
449
450
451
		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);

452
		if ( config_gpu ) {
453
454
			if ( gctx == NULL ) {
				gctx = setup_gpu(config_nosfac, &image,
Thomas White's avatar
Thomas White committed
455
				                 intensities, gpu_dev);
456
			}
457
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
458
		} else {
Thomas White's avatar
Thomas White committed
459
460
			get_diffraction(&image, na, nb, nc, intensities, phases,
			                cell, !config_nowater, grad);
461
		}
462
		if ( image.data == NULL ) {
463
464
465
466
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

467
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
468

Thomas White's avatar
Thomas White committed
469
		if ( config_nearbragg ) {
470
			find_projected_peaks(&image, cell, 0, 0.1);
471
472
			output_intensities(&image, cell, NULL, 0, 1, 0, stdout,
			                   0, 0.1);
Thomas White's avatar
Thomas White committed
473
			free(image.cpeaks);
Thomas White's avatar
Thomas White committed
474
475
		}

476
		if ( powder_fn != NULL ) {
477
478
479
480
481
482
483

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
484
				powder[x+w*y] += (double)image.data[x+w*y];
485
486
487
488
			}
			}

			if ( !(ndone % 10) ) {
489
				hdf5_write(powder_fn, powder,
490
				           image.width, image.height,
491
				           H5T_NATIVE_DOUBLE);
492
493
494
			}
		}

495
496
497
498
		if ( !config_noimages ) {

			char filename[1024];

499
500
501
502
503
504
505
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

506
			number++;
Thomas White's avatar
Thomas White committed
507

508
509
			/* Write the output file */
			hdf5_write(filename, image.data,
510
			           image.width, image.height, H5T_NATIVE_FLOAT);
511
512

		}
513

514
515
516
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
517
		cell_free(cell);
Thomas White's avatar
Thomas White committed
518

519
skip:
Thomas White's avatar
Thomas White committed
520
		ndone++;
521

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

524
	} while ( !done );
525

526
527
528
529
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
530
531
	free(image.det->panels);
	free(image.det);
532
	free(image.beam);
533
	free(powder);
Thomas White's avatar
Thomas White committed
534
	cell_free(input_cell);
535
	free(intensities);
Thomas White's avatar
Thomas White committed
536
537
	free(outfile);
	free(filename);
538

539
	return 0;
540
}