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


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


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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	}

276
277
278
279
280
281
282
283
284
285
	/* 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);
	}

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

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

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

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

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

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

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

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

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

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

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

397
398
	do {

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

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

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

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

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

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

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

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

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

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

			int x, y, w;

			w = image.width;

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

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

475
476
477
478
		if ( !config_noimages ) {

			char filename[1024];

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

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

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

		}
493

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

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

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

503
	} while ( !done );
504

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

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

518
	return 0;
519
}