pattern_sim.c 13.8 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"
Thomas White's avatar
Thomas White committed
83
84
85
86
87
"\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"
88
"     --no-water            Do not simulate water background.\n"
Thomas White's avatar
Thomas White committed
89
"     --no-noise            Do not calculate Poisson noise.\n"
90
);
Thomas White's avatar
Thomas White committed
91
92
93
}


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


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

	} while ( 1 );
169
170
171
172
173
}


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

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

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

Thomas White's avatar
Thomas White committed
227
		switch (c) {
Thomas White's avatar
Thomas White committed
228
		case 'h' :
Thomas White's avatar
Thomas White committed
229
230
			show_help(argv[0]);
			return 0;
231

Thomas White's avatar
Thomas White committed
232
		case 'r' :
233
234
235
			config_randomquat = 1;
			break;

Thomas White's avatar
Thomas White committed
236
		case 'n' :
237
238
239
			n_images = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
240
		case 'i' :
241
242
243
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
244
		case 't' :
245
246
247
			grad_str = strdup(optarg);
			break;

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

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

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

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

264
265
266
267
		case 'b' :
			beamfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
268
		case 0 :
Thomas White's avatar
Thomas White committed
269
			break;
270

Thomas White's avatar
Thomas White committed
271
		default :
Thomas White's avatar
Thomas White committed
272
273
			return 1;
		}
274
275
276

	}

277
278
279
280
281
282
283
284
285
286
	/* FIXME: Make this optional */
	if ( 1 ) {
		FILE *fh;
		unsigned int seed;
		fh = fopen("/dev/urandom", "r");
		fread(&seed, sizeof(seed), 1, fh);
		fclose(fh);
		srand(seed);
	}

287
288
289
290
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

291
292
293
294
295
296
297
298
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

Thomas White's avatar
Thomas White committed
299
	if ( config_simdetails ) {
300
		show_details();
Thomas White's avatar
Thomas White committed
301
302
303
		return 0;
	}

304
305
306
307
308
309
	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;
	}

310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
	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
332
333
334
335
336
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

337
338
339
340
341
342
	if ( beamfile == NULL ) {
		ERROR("You need to specify a beam parameter file"
		      " with --beam\n");
		return 1;
	}

343
344
345
346
347
348
349
350
	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
351
		ReflItemList *items;
352
353
354
355
356
		if ( grad == GRADIENT_PHASED ) {
			phases = new_list_phase();
		} else {
			phases = NULL;
		}
Thomas White's avatar
Thomas White committed
357
358
		intensities = new_list_intensity();
		phases = new_list_phase();
359
360
		items = read_reflections(intfile, intensities, phases,
		                         NULL, NULL);
361
		free(intfile);
Thomas White's avatar
Thomas White committed
362
		delete_items(items);
363
364
	}

365
366
367
368
369
370
371
	image.det = get_detector_geometry(geometry);
	if ( image.det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

372
373
374
375
376
377
378
	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
379
	/* Define image parameters */
380
381
	image.width = image.det->max_x + 1;
	image.height = image.det->max_y + 1;
382
	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
383
	cell = load_cell_from_pdb(filename);
Thomas White's avatar
Thomas White committed
384
385
386
	if ( cell == NULL ) {
		exit(1);
	}
387
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
388
389
	image.features = NULL;
	image.flags = NULL;
390
391
	image.f0 = 1.0;
	image.f0_available = 1;
Thomas White's avatar
Thomas White committed
392

393
394
	powder = calloc(image.width*image.height, sizeof(*powder));

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

398
399
	do {

Thomas White's avatar
Thomas White committed
400
		int na, nb, nc;
401
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
402

403
404
405
		//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
406
407
408
		na = 24;
		nb = 24;
		nc = 40;
Thomas White's avatar
Thomas White committed
409

410
411
412
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
			image.orientation = random_quaternion();
413
		} else {
414
			image.orientation = read_quaternion();
415
416
		}

Thomas White's avatar
Thomas White committed
417
418
419
420
		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
421
		if ( !quaternion_valid(image.orientation) ) {
422
			ERROR("Orientation modulus is not zero!\n");
423
424
425
426
427
428
			return 1;
		}

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

430
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
431
432
433
		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);

434
		if ( config_gpu ) {
435
436
			if ( gctx == NULL ) {
				gctx = setup_gpu(config_nosfac, &image,
Thomas White's avatar
Thomas White committed
437
				                 intensities);
438
			}
439
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
440
		} else {
Thomas White's avatar
Thomas White committed
441
442
			get_diffraction(&image, na, nb, nc, intensities, phases,
			                cell, !config_nowater, grad);
443
		}
444
		if ( image.data == NULL ) {
445
446
447
448
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

449
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
450

Thomas White's avatar
Thomas White committed
451
		if ( config_nearbragg ) {
452
			find_projected_peaks(&image, cell, 0, 0.1);
453
454
			output_intensities(&image, cell, NULL, 0, 1, 0, stdout,
			                   0, 0.1);
Thomas White's avatar
Thomas White committed
455
456
		}

457
		if ( powder_fn != NULL ) {
458
459
460
461
462
463
464

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
465
				powder[x+w*y] += (double)image.data[x+w*y];
466
467
468
469
			}
			}

			if ( !(ndone % 10) ) {
470
				hdf5_write(powder_fn, powder,
471
				           image.width, image.height,
472
				           H5T_NATIVE_DOUBLE);
473
474
475
			}
		}

476
477
478
479
		if ( !config_noimages ) {

			char filename[1024];

480
481
482
483
484
485
486
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

487
			number++;
Thomas White's avatar
Thomas White committed
488

489
490
			/* Write the output file */
			hdf5_write(filename, image.data,
491
			           image.width, image.height, H5T_NATIVE_FLOAT);
492
493

		}
494

495
496
497
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
498

499
skip:
Thomas White's avatar
Thomas White committed
500
		ndone++;
501

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

504
	} while ( !done );
505

506
507
508
509
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
510
511
	free(image.det->panels);
	free(image.det);
512
	free(image.beam);
513
	free(powder);
514
	cell_free(cell);
515
	free(intensities);
Thomas White's avatar
Thomas White committed
516
517
	free(outfile);
	free(filename);
518

519
	return 0;
520
}