pattern_sim.c 15.7 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-2011 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"
32
#include "beam-parameters.h"
Thomas White's avatar
Thomas White committed
33
#include "symmetry.h"
Thomas White's avatar
Thomas White committed
34
35
#include "reflist.h"
#include "reflist-utils.h"
36
37


Thomas White's avatar
Thomas White committed
38
static void show_help(const char *s)
39
{
40
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
41
	printf(
Thomas White's avatar
Thomas White committed
42
"Simulate diffraction patterns from small crystals probed with femtosecond\n"
Thomas White's avatar
Thomas White committed
43
44
"pulses of X-rays from a free electron laser.\n"
"\n"
45
" -h, --help                Display this help message.\n"
Thomas White's avatar
Thomas White committed
46
"\n"
47
48
49
" -p, --pdb=<file>          PDB file from which to get the unit cell.\n"
"                            (The actual Bragg intensities come from the\n"
"                            intensities file)\n"
50
"     --simulation-details  Show technical details of the simulation.\n"
51
"     --gpu                 Use the GPU to speed up the calculation.\n"
Thomas White's avatar
Thomas White committed
52
53
"     --gpu-dev=<n>         Use GPU device <n>.  Omit this option to see the\n"
"                            available devices.\n"
54
55
" -g, --geometry=<file>     Get detector geometry from file.\n"
" -b, --beam=<file>         Get beam parameters from file.\n"
56
57
58
"\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"
71
" -y, --symmetry=<sym>      The symmetry of the intensities file.\n"
Thomas White's avatar
Thomas White committed
72
" -t, --gradients=<method>  Use <method> for the calculation of shape\n"
73
74
"                            transform intensities.  Choose from:\n"
"                             mosaic      : Take the intensity of the nearest\n"
Thomas White's avatar
Thomas White committed
75
76
77
"                                           Bragg position.  This is the\n"
"                                           fastest method and the only one\n"
"                                           supported on the GPU.\n"
78
79
"                             interpolate : Interpolate trilinearly between\n"
"                                           six adjacent Bragg intensities.\n"
Thomas White's avatar
Thomas White committed
80
81
"                                           This method has intermediate\n"
"                                           accuracy.\n"
82
"                             phased      : As 'interpolate', but take phase\n"
Thomas White's avatar
Thomas White committed
83
84
85
"                                           values into account.  This is the\n"
"                                           most accurate method, but the\n"
"                                           slowest.\n"
86
87
88
"     --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"
89
90
91
92
93
"     --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
94
95
96
97
98
"\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"
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
125
126
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
127
"I(q) = I0 * r^2 * I_latt(q) * I_mol(q) * S\n"
Thomas White's avatar
Thomas White committed
128
129
130
131
"\n"
"I0 = number of photons per unit area in the incident beam\n"
" r = Thomson radius\n"
" S = solid angle of corresponding pixel\n"
132
133
134
"\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
135
136
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
137
138
139
"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"
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
168
169
170
}


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

	} while ( 1 );
175
176
177
}


178
179
180
181
182
183
184
185
186
187
188
189
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;
}


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

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

251
	/* Short options */
252
	while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:g:b:y:",
Thomas White's avatar
Thomas White committed
253
	                        longopts, NULL)) != -1) {
254

Thomas White's avatar
Thomas White committed
255
		switch (c) {
Thomas White's avatar
Thomas White committed
256
		case 'h' :
Thomas White's avatar
Thomas White committed
257
258
			show_help(argv[0]);
			return 0;
259

Thomas White's avatar
Thomas White committed
260
		case 'r' :
261
262
263
			config_randomquat = 1;
			break;

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

Thomas White's avatar
Thomas White committed
272
		case 'i' :
273
274
275
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
276
		case 't' :
277
278
279
			grad_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
280
		case 'p' :
281
282
283
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
284
		case 'o' :
285
286
287
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
288
		case 'w' :
289
290
291
			powder_fn = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
292
		case 'g' :
Thomas White's avatar
Thomas White committed
293
294
295
			geometry = strdup(optarg);
			break;

296
297
298
299
		case 'b' :
			beamfile = strdup(optarg);
			break;

300
301
302
303
		case 'y' :
			sym = strdup(optarg);
			break;

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

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

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

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

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

	}

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

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

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

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

361
362
	if ( sym == NULL ) sym = strdup("1");

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

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

395
396
397
398
399
400
	if ( beamfile == NULL ) {
		ERROR("You need to specify a beam parameter file"
		      " with --beam\n");
		return 1;
	}

401
	if ( intfile == NULL ) {
402

403
404
405
406
407
408
		/* 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;
409
410
		flags = NULL;

411
	} else {
412

Thomas White's avatar
Thomas White committed
413
414
415
416
		RefList *reflections;

		reflections = read_reflections(intfile);
		free(intfile);
417

418
		if ( grad == GRADIENT_PHASED ) {
Thomas White's avatar
Thomas White committed
419
			phases = phases_from_list(reflections);
420
421
422
		} else {
			phases = NULL;
		}
Thomas White's avatar
Thomas White committed
423
424
425
		intensities = intensities_from_list(reflections);
		phases = phases_from_list(reflections);
		flags = flags_from_list(reflections);
426

Thomas White's avatar
Thomas White committed
427
		/* Check that the intensities have the correct symmetry */
Thomas White's avatar
Thomas White committed
428
		if ( check_list_symmetry(reflections, sym) ) {
Thomas White's avatar
Thomas White committed
429
430
431
432
433
			ERROR("The input reflection list does not appear to"
			      " have symmetry %s\n", sym);
			return 1;
		}

Thomas White's avatar
Thomas White committed
434
435
		reflist_free(reflections);

436
437
	}

438
439
440
441
442
443
444
	image.det = get_detector_geometry(geometry);
	if ( image.det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

445
446
447
448
449
450
451
	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
452
	/* Define image parameters */
Thomas White's avatar
Thomas White committed
453
454
	image.width = image.det->max_fs + 1;
	image.height = image.det->max_ss + 1;
455
	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
Thomas White's avatar
Thomas White committed
456
457
458
459

	/* Load unit cell */
	input_cell = load_cell_from_pdb(filename);
	if ( input_cell == NULL ) {
Thomas White's avatar
Thomas White committed
460
461
		exit(1);
	}
Thomas White's avatar
Thomas White committed
462
463

	/* Initialise stuff */
464
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
465
466
	image.features = NULL;
	image.flags = NULL;
Thomas White's avatar
Thomas White committed
467
468
	image.i0 = 1.0;
	image.i0_available = 1;
Thomas White's avatar
Thomas White committed
469

470
471
	powder = calloc(image.width*image.height, sizeof(*powder));

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

475
476
	do {

Thomas White's avatar
Thomas White committed
477
		int na, nb, nc;
478
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
479
		UnitCell *cell;
Thomas White's avatar
Thomas White committed
480

481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
		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
499

500
501
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
Thomas White's avatar
Thomas White committed
502
			orientation = random_quaternion();
503
		} else {
Thomas White's avatar
Thomas White committed
504
			orientation = read_quaternion();
505
506
		}

Thomas White's avatar
Thomas White committed
507
		STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
Thomas White's avatar
Thomas White committed
508
509
		       orientation.w, orientation.x,
		       orientation.y, orientation.z);
Thomas White's avatar
Thomas White committed
510

Thomas White's avatar
Thomas White committed
511
		if ( !quaternion_valid(orientation) ) {
512
			ERROR("Orientation modulus is not zero!\n");
513
514
515
			return 1;
		}

Thomas White's avatar
Thomas White committed
516
517
		cell = cell_rotate(input_cell, orientation);

518
519
520
		/* Ensure no residual information */
		image.data = NULL;
		image.twotheta = NULL;
521

522
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
Thomas White's avatar
Thomas White committed
523
524
		STATUS("Particle size = %i x %i x %i"
		       " ( = %5.2f x %5.2f x %5.2f nm)\n",
525
526
	               na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9);

527
		if ( config_gpu ) {
528
			if ( gctx == NULL ) {
529
				gctx = setup_gpu(config_nosfac,
530
531
				                 intensities, flags, sym,
				                 gpu_dev);
532
			}
533
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
534
		} else {
Thomas White's avatar
Thomas White committed
535
			get_diffraction(&image, na, nb, nc, intensities, phases,
Thomas White's avatar
Thomas White committed
536
			                flags, cell, grad, sym);
537
		}
538
		if ( image.data == NULL ) {
539
540
541
542
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

543
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
544

545
		if ( powder_fn != NULL ) {
546
547
548
549
550
551
552

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
553
				powder[x+w*y] += (double)image.data[x+w*y];
554
555
556
557
			}
			}

			if ( !(ndone % 10) ) {
558
				hdf5_write(powder_fn, powder,
559
				           image.width, image.height,
560
				           H5T_NATIVE_DOUBLE);
561
562
563
			}
		}

564
565
566
567
		if ( !config_noimages ) {

			char filename[1024];

568
569
570
571
572
573
574
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

575
			number++;
Thomas White's avatar
Thomas White committed
576

577
578
			/* Write the output file */
			hdf5_write(filename, image.data,
579
			           image.width, image.height, H5T_NATIVE_FLOAT);
580
581

		}
582

583
584
585
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
586
		cell_free(cell);
Thomas White's avatar
Thomas White committed
587

588
skip:
Thomas White's avatar
Thomas White committed
589
		ndone++;
590

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

593
	} while ( !done );
594

595
596
597
598
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
599
600
	free(image.det->panels);
	free(image.det);
601
	free(image.beam);
602
	free(powder);
Thomas White's avatar
Thomas White committed
603
	cell_free(input_cell);
604
	free(intensities);
Thomas White's avatar
Thomas White committed
605
606
	free(outfile);
	free(filename);
607
	free(sym);
608

609
	return 0;
610
}