pattern_sim.c 14.5 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"
Thomas White's avatar
Thomas White committed
88
89
90
91
92
"\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"
93
"     --no-water            Do not simulate water background.\n"
Thomas White's avatar
Thomas White committed
94
"     --no-noise            Do not calculate Poisson noise.\n"
95
);
Thomas White's avatar
Thomas White committed
96
97
98
}


99
static void show_details()
Thomas White's avatar
Thomas White committed
100
{
Thomas White's avatar
Thomas White committed
101
102
103
104
105
	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"
106
107
"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
108
"\n"
109
110
111
"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
112
113
114
115
"\n"
"na = number of unit cells in 'a' direction (likewise nb, nc)\n"
" q = reciprocal vector (1/d convention, not 2pi/d)\n"
"\n"
116
117
118
"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
119
"\n"
120
121
"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"
122
"coherent, elastic part of the diffuse scattering from the water jet only.\n"
Thomas White's avatar
Thomas White committed
123
124
125
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
126
"I(q) = I0 * r^2 * I_latt(q) * I_mol(q) * S\n"
Thomas White's avatar
Thomas White committed
127
128
129
130
"\n"
"I0 = number of photons per unit area in the incident beam\n"
" r = Thomson radius\n"
" S = solid angle of corresponding pixel\n"
131
132
133
"\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
134
135
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
136
137
138
"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"
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
168
169
}


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

	} while ( 1 );
174
175
176
177
178
}


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

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

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

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

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

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

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

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

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

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

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

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

274
275
276
277
		case 'b' :
			beamfile = strdup(optarg);
			break;

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

Thomas White's avatar
Thomas White committed
282
		case 0 :
Thomas White's avatar
Thomas White committed
283
			break;
284

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

	}

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

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

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

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

317
318
319
320
321
322
	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;
	}

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

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

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

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

385
386
387
388
389
390
391
	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
392
	/* Define image parameters */
393
394
	image.width = image.det->max_x + 1;
	image.height = image.det->max_y + 1;
395
	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
Thomas White's avatar
Thomas White committed
396
397
398
399

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

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

410
411
	powder = calloc(image.width*image.height, sizeof(*powder));

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

415
416
	do {

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

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

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

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

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

Thomas White's avatar
Thomas White committed
444
445
		cell = cell_rotate(input_cell, orientation);

446
447
448
		/* Ensure no residual information */
		image.data = NULL;
		image.twotheta = NULL;
449

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

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

469
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
470

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

478
		if ( powder_fn != NULL ) {
479
480
481
482
483
484
485

			int x, y, w;

			w = image.width;

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

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

497
498
499
500
		if ( !config_noimages ) {

			char filename[1024];

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

508
			number++;
Thomas White's avatar
Thomas White committed
509

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

		}
515

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

521
skip:
Thomas White's avatar
Thomas White committed
522
		ndone++;
523

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

526
	} while ( !done );
527

528
529
530
531
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

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

541
	return 0;
542
}