pattern_sim.c 15.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 "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"
Thomas White's avatar
Thomas White committed
51
52
"     --gpu-dev=<n>         Use GPU device <n>.  Omit this option to see the\n"
"                            available devices.\n"
53
54
" -g, --geometry=<file>     Get detector geometry from file.\n"
" -b, --beam=<file>         Get beam parameters from file.\n"
55
56
57
58
"\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"
59
60
61
62
63
64
" -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
65
" -r, --random-orientation  Use a randomly generated orientation\n"
66
"                            (a new orientation will be used for each image).\n"
67
68
"     --powder=<file>       Write a summed pattern of all images simulated by\n"
"                            this invocation as the given filename.\n"
69
" -i, --intensities=<file>  Specify file containing reflection intensities\n"
Thomas White's avatar
Thomas White committed
70
"                            (and phases) to use.\n"
Thomas White's avatar
Thomas White committed
71
" -t, --gradients=<method>  Use <method> for the calculation of shape\n"
72
73
"                            transform intensities.  Choose from:\n"
"                             mosaic      : Take the intensity of the nearest\n"
Thomas White's avatar
Thomas White committed
74
75
76
"                                           Bragg position.  This is the\n"
"                                           fastest method and the only one\n"
"                                           supported on the GPU.\n"
77
78
"                             interpolate : Interpolate trilinearly between\n"
"                                           six adjacent Bragg intensities.\n"
Thomas White's avatar
Thomas White committed
79
80
"                                           This method has intermediate\n"
"                                           accuracy.\n"
81
"                             phased      : As 'interpolate', but take phase\n"
Thomas White's avatar
Thomas White committed
82
83
84
"                                           values into account.  This is the\n"
"                                           most accurate method, but the\n"
"                                           slowest.\n"
85
86
87
"     --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"
88
89
90
91
92
"     --min-size=<s>        Generate random crystal sizes using <s> as the\n"
"                            minimum crystal size in nm. --max-size is also\n"
"                            required.\n"
"     --max-size=<s>        Use <s> as the maximum crystal size in nm.\n"
"                            --min-size is also required.\n"
Thomas White's avatar
Thomas White committed
93
94
95
96
97
"\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"
98
"     --no-water            Do not simulate water background.\n"
Thomas White's avatar
Thomas White committed
99
"     --no-noise            Do not calculate Poisson noise.\n"
100
);
Thomas White's avatar
Thomas White committed
101
102
103
}


104
static void show_details()
Thomas White's avatar
Thomas White committed
105
{
Thomas White's avatar
Thomas White committed
106
107
108
109
110
	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"
111
112
"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
113
"\n"
114
115
116
"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
117
118
119
120
"\n"
"na = number of unit cells in 'a' direction (likewise nb, nc)\n"
" q = reciprocal vector (1/d convention, not 2pi/d)\n"
"\n"
121
122
123
"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
124
"\n"
125
126
"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"
127
"coherent, elastic part of the diffuse scattering from the water jet only.\n"
Thomas White's avatar
Thomas White committed
128
129
130
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
131
"I(q) = I0 * r^2 * I_latt(q) * I_mol(q) * S\n"
Thomas White's avatar
Thomas White committed
132
133
134
135
"\n"
"I0 = number of photons per unit area in the incident beam\n"
" r = Thomson radius\n"
" S = solid angle of corresponding pixel\n"
136
137
138
"\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
139
140
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
141
142
143
"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"
144
);
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
}


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 {
175
			ERROR("Invalid rotation '%s'\n", line);
176
177
178
		}

	} while ( 1 );
179
180
181
}


182
183
184
185
186
187
188
189
190
191
192
193
static int random_ncells(double len, double min_nm, double max_nm)
{
	int min_cells, max_cells, wid;

	min_cells = min_nm*1e-9 / len;
	max_cells = max_nm*1e-9 / len;
	wid = max_cells - min_cells;

	return wid*random()/RAND_MAX + min_cells;
}


194
195
int main(int argc, char *argv[])
{
196
	int c;
197
	struct image image;
198
	struct gpu_context *gctx = NULL;
199
	double *powder;
200
201
	char *intfile = NULL;
	double *intensities;
202
	char *rval;
203
	double *phases;
Thomas White's avatar
Thomas White committed
204
	int config_simdetails = 0;
205
206
	int config_nearbragg = 0;
	int config_randomquat = 0;
207
	int config_noimages = 0;
208
	int config_nowater = 0;
Thomas White's avatar
Thomas White committed
209
	int config_nonoise = 0;
210
	int config_nosfac = 0;
211
	int config_gpu = 0;
212
	int config_random = 0;
213
	char *powder_fn = NULL;
214
	char *filename = NULL;
215
	char *grad_str = NULL;
216
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
217
	char *geometry = NULL;
218
	char *beamfile = NULL;
219
	GradientMethod grad;
Thomas White's avatar
Thomas White committed
220
221
222
	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 */
223
	int done = 0;
Thomas White's avatar
Thomas White committed
224
225
	UnitCell *input_cell;
	struct quaternion orientation;
Thomas White's avatar
Thomas White committed
226
	int gpu_dev = -1;
227
228
229
	int random_size = 0;
	double min_size = 0.0;
	double max_size = 0.0;
Thomas White's avatar
Thomas White committed
230

231
	/* Long options */
Thomas White's avatar
Thomas White committed
232
	const struct option longopts[] = {
233
234
		{"help",               0, NULL,               'h'},
		{"simulation-details", 0, &config_simdetails,  1},
235
		{"gpu",                0, &config_gpu,         1},
236
237
238
		{"near-bragg",         0, &config_nearbragg,   1},
		{"random-orientation", 0, NULL,               'r'},
		{"number",             1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
239
		{"no-images",          0, &config_noimages,    1},
240
		{"no-water",           0, &config_nowater,     1},
Thomas White's avatar
Thomas White committed
241
		{"no-noise",           0, &config_nonoise,     1},
242
		{"intensities",        1, NULL,               'i'},
243
		{"powder",             1, NULL,               'w'},
Thomas White's avatar
Thomas White committed
244
		{"gradients",          1, NULL,               't'},
245
		{"pdb",                1, NULL,               'p'},
246
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
247
		{"geometry",           1, NULL,               'g'},
248
		{"beam",               1, NULL,               'b'},
249
		{"really-random",      0, &config_random,      1},
Thomas White's avatar
Thomas White committed
250
		{"gpu-dev",            1, NULL,                2},
251
252
		{"min-size",           1, NULL,                3},
		{"max-size",           1, NULL,                4},
253
		{0, 0, NULL, 0}
Thomas White's avatar
Thomas White committed
254
	};
255

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

Thomas White's avatar
Thomas White committed
260
		switch (c) {
Thomas White's avatar
Thomas White committed
261
		case 'h' :
Thomas White's avatar
Thomas White committed
262
263
			show_help(argv[0]);
			return 0;
264

Thomas White's avatar
Thomas White committed
265
		case 'r' :
266
267
268
			config_randomquat = 1;
			break;

Thomas White's avatar
Thomas White committed
269
		case 'n' :
270
271
272
273
274
			n_images = strtol(optarg, &rval, 10);
			if ( *rval != '\0' ) {
				ERROR("Invalid number of images.\n");
				return 1;
			}
275
276
			break;

Thomas White's avatar
Thomas White committed
277
		case 'i' :
278
279
280
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
281
		case 't' :
282
283
284
			grad_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
285
		case 'p' :
286
287
288
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
289
		case 'o' :
290
291
292
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
293
		case 'w' :
294
295
296
			powder_fn = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
297
		case 'g' :
Thomas White's avatar
Thomas White committed
298
299
300
			geometry = strdup(optarg);
			break;

301
302
303
304
		case 'b' :
			beamfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
305
306
307
308
		case 2 :
			gpu_dev = atoi(optarg);
			break;

309
		case 3 :
310
311
312
313
314
			min_size = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid minimum size.\n");
				return 1;
			}
315
316
317
318
			random_size++;
			break;

		case 4 :
319
320
321
322
323
			max_size = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid maximum size.\n");
				return 1;
			}
324
325
326
			random_size++;
			break;

Thomas White's avatar
Thomas White committed
327
		case 0 :
Thomas White's avatar
Thomas White committed
328
			break;
329

Thomas White's avatar
Thomas White committed
330
		default :
Thomas White's avatar
Thomas White committed
331
332
			return 1;
		}
333
334
335

	}

336
	if ( config_random ) {
337
338
339
340
341
342
343
344
		FILE *fh;
		unsigned int seed;
		fh = fopen("/dev/urandom", "r");
		fread(&seed, sizeof(seed), 1, fh);
		fclose(fh);
		srand(seed);
	}

345
346
347
348
349
	if ( random_size == 1 ) {
		ERROR("You must specify both --min-size and --max-size.\n");
		return 1;
	}

350
351
352
353
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

354
355
356
357
358
359
360
361
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

Thomas White's avatar
Thomas White committed
362
	if ( config_simdetails ) {
363
		show_details();
Thomas White's avatar
Thomas White committed
364
365
366
		return 0;
	}

367
368
369
370
371
372
	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;
	}

373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
	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
395
396
397
398
399
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

400
401
402
403
404
405
	if ( beamfile == NULL ) {
		ERROR("You need to specify a beam parameter file"
		      " with --beam\n");
		return 1;
	}

406
407
408
409
410
411
412
413
	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
414
		ReflItemList *items;
415
416
417
418
419
		if ( grad == GRADIENT_PHASED ) {
			phases = new_list_phase();
		} else {
			phases = NULL;
		}
Thomas White's avatar
Thomas White committed
420
421
		intensities = new_list_intensity();
		phases = new_list_phase();
422
423
		items = read_reflections(intfile, intensities, phases,
		                         NULL, NULL);
424
		free(intfile);
Thomas White's avatar
Thomas White committed
425
		delete_items(items);
426
427
	}

428
429
430
431
432
433
434
	image.det = get_detector_geometry(geometry);
	if ( image.det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

435
436
437
438
439
440
441
	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
442
	/* Define image parameters */
443
444
	image.width = image.det->max_x + 1;
	image.height = image.det->max_y + 1;
445
	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
Thomas White's avatar
Thomas White committed
446
447
448
449

	/* Load unit cell */
	input_cell = load_cell_from_pdb(filename);
	if ( input_cell == NULL ) {
Thomas White's avatar
Thomas White committed
450
451
		exit(1);
	}
Thomas White's avatar
Thomas White committed
452
453

	/* Initialise stuff */
454
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
455
456
	image.features = NULL;
	image.flags = NULL;
457
458
	image.f0 = 1.0;
	image.f0_available = 1;
Thomas White's avatar
Thomas White committed
459

460
461
	powder = calloc(image.width*image.height, sizeof(*powder));

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

465
466
	do {

Thomas White's avatar
Thomas White committed
467
		int na, nb, nc;
468
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
469
		UnitCell *cell;
Thomas White's avatar
Thomas White committed
470

471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
		if ( random_size ) {

			double alen, blen, clen, dis;

			cell_get_parameters(input_cell, &alen, &blen, &clen,
			                    &dis, &dis, &dis);

			na = random_ncells(alen, min_size, max_size);
			nb = random_ncells(blen, min_size, max_size);
			nc = random_ncells(clen, min_size, max_size);

		} else {

			na = 8;
			nb = 8;
			nc = 8;

		}
Thomas White's avatar
Thomas White committed
489

490
491
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
Thomas White's avatar
Thomas White committed
492
			orientation = random_quaternion();
493
		} else {
Thomas White's avatar
Thomas White committed
494
			orientation = read_quaternion();
495
496
		}

Thomas White's avatar
Thomas White committed
497
		STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
Thomas White's avatar
Thomas White committed
498
499
		       orientation.w, orientation.x,
		       orientation.y, orientation.z);
Thomas White's avatar
Thomas White committed
500

Thomas White's avatar
Thomas White committed
501
		if ( !quaternion_valid(orientation) ) {
502
			ERROR("Orientation modulus is not zero!\n");
503
504
505
			return 1;
		}

Thomas White's avatar
Thomas White committed
506
507
		cell = cell_rotate(input_cell, orientation);

508
509
510
		/* Ensure no residual information */
		image.data = NULL;
		image.twotheta = NULL;
511

512
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
513
514
515
		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);

516
		if ( config_gpu ) {
517
518
			if ( gctx == NULL ) {
				gctx = setup_gpu(config_nosfac, &image,
Thomas White's avatar
Thomas White committed
519
				                 intensities, gpu_dev);
520
			}
521
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
522
		} else {
Thomas White's avatar
Thomas White committed
523
524
			get_diffraction(&image, na, nb, nc, intensities, phases,
			                cell, !config_nowater, grad);
525
		}
526
		if ( image.data == NULL ) {
527
528
529
530
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

531
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
532

Thomas White's avatar
Thomas White committed
533
		if ( config_nearbragg ) {
534
			find_projected_peaks(&image, cell, 0, 0.1);
535
536
			output_intensities(&image, cell, NULL, 0, 1, 0, stdout,
			                   0, 0.1);
Thomas White's avatar
Thomas White committed
537
			free(image.cpeaks);
Thomas White's avatar
Thomas White committed
538
539
		}

540
		if ( powder_fn != NULL ) {
541
542
543
544
545
546
547

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
548
				powder[x+w*y] += (double)image.data[x+w*y];
549
550
551
552
			}
			}

			if ( !(ndone % 10) ) {
553
				hdf5_write(powder_fn, powder,
554
				           image.width, image.height,
555
				           H5T_NATIVE_DOUBLE);
556
557
558
			}
		}

559
560
561
562
		if ( !config_noimages ) {

			char filename[1024];

563
564
565
566
567
568
569
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

570
			number++;
Thomas White's avatar
Thomas White committed
571

572
573
			/* Write the output file */
			hdf5_write(filename, image.data,
574
			           image.width, image.height, H5T_NATIVE_FLOAT);
575
576

		}
577

578
579
580
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
581
		cell_free(cell);
Thomas White's avatar
Thomas White committed
582

583
skip:
Thomas White's avatar
Thomas White committed
584
		ndone++;
585

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

588
	} while ( !done );
589

590
591
592
593
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
594
595
	free(image.det->panels);
	free(image.det);
596
	free(image.beam);
597
	free(powder);
Thomas White's avatar
Thomas White committed
598
	cell_free(input_cell);
599
	free(intensities);
Thomas White's avatar
Thomas White committed
600
601
	free(outfile);
	free(filename);
602

603
	return 0;
604
}