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


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


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

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


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

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

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

Thomas White's avatar
Thomas White committed
223
		switch (c) {
Thomas White's avatar
Thomas White committed
224
		case 'h' :
Thomas White's avatar
Thomas White committed
225
226
			show_help(argv[0]);
			return 0;
227

Thomas White's avatar
Thomas White committed
228
		case 'r' :
229
230
231
			config_randomquat = 1;
			break;

Thomas White's avatar
Thomas White committed
232
		case 'n' :
233
234
235
			n_images = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
236
		case 'i' :
237
238
239
			intfile = strdup(optarg);
			break;

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

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

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

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

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

Thomas White's avatar
Thomas White committed
260
		case 0 :
Thomas White's avatar
Thomas White committed
261
			break;
262

Thomas White's avatar
Thomas White committed
263
		default :
Thomas White's avatar
Thomas White committed
264
265
			return 1;
		}
266
267
268

	}

269
270
271
272
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

273
274
275
276
277
278
279
280
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

Thomas White's avatar
Thomas White committed
281
	if ( config_simdetails ) {
282
		show_details();
Thomas White's avatar
Thomas White committed
283
284
285
		return 0;
	}

286
287
288
289
290
291
	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;
	}

292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
	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
314
315
316
317
318
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

319
320
321
322
323
324
325
326
	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
327
		ReflItemList *items;
328
329
330
331
332
		if ( grad == GRADIENT_PHASED ) {
			phases = new_list_phase();
		} else {
			phases = NULL;
		}
Thomas White's avatar
Thomas White committed
333
334
		intensities = new_list_intensity();
		phases = new_list_phase();
335
336
		items = read_reflections(intfile, intensities, phases,
		                         NULL, NULL);
337
		free(intfile);
Thomas White's avatar
Thomas White committed
338
		delete_items(items);
339
340
	}

Thomas White's avatar
Thomas White committed
341
	/* Define image parameters */
Thomas White's avatar
Thomas White committed
342
343
	image.width = 1024;
	image.height = 1024;
344
	image.lambda = ph_en_to_lambda(eV_to_J(PHOTON_ENERGY)); /* Wavelength */
345
	cell = load_cell_from_pdb(filename);
Thomas White's avatar
Thomas White committed
346
347
348
	if ( cell == NULL ) {
		exit(1);
	}
349
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
350
351
	image.features = NULL;
	image.flags = NULL;
352
353
	image.f0 = 1.0;
	image.f0_available = 1;
Thomas White's avatar
Thomas White committed
354

Thomas White's avatar
Thomas White committed
355
356
357
358
359
360
	image.det = get_detector_geometry(geometry);
	if ( image.det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);
361

362
363
	powder = calloc(image.width*image.height, sizeof(*powder));

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

367
368
	do {

Thomas White's avatar
Thomas White committed
369
		int na, nb, nc;
370
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
371

372
373
374
		//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
375
376
377
		na = 24;
		nb = 24;
		nc = 40;
Thomas White's avatar
Thomas White committed
378

379
380
381
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
			image.orientation = random_quaternion();
382
		} else {
383
			image.orientation = read_quaternion();
384
385
		}

Thomas White's avatar
Thomas White committed
386
387
388
389
		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
390
		if ( !quaternion_valid(image.orientation) ) {
391
			ERROR("Orientation modulus is not zero!\n");
392
393
394
395
396
397
			return 1;
		}

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

399
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
400
401
402
		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);

403
		if ( config_gpu ) {
404
405
			if ( gctx == NULL ) {
				gctx = setup_gpu(config_nosfac, &image,
Thomas White's avatar
Thomas White committed
406
				                 intensities);
407
			}
408
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
409
		} else {
Thomas White's avatar
Thomas White committed
410
411
			get_diffraction(&image, na, nb, nc, intensities, phases,
			                cell, !config_nowater, grad);
412
		}
413
		if ( image.data == NULL ) {
414
415
416
417
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

418
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
419

Thomas White's avatar
Thomas White committed
420
		if ( config_nearbragg ) {
421
422
			find_projected_peaks(&image, cell, 0, 0.1);
			output_intensities(&image, cell, NULL, 0, 1, 0, 0, 0.1);
Thomas White's avatar
Thomas White committed
423
424
		}

425
		if ( powder_fn != NULL ) {
426
427
428
429
430
431
432

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
433
				powder[x+w*y] += (double)image.data[x+w*y];
434
435
436
437
			}
			}

			if ( !(ndone % 10) ) {
438
				hdf5_write(powder_fn, powder,
439
				           image.width, image.height,
440
				           H5T_NATIVE_DOUBLE);
441
442
443
			}
		}

444
445
446
447
		if ( !config_noimages ) {

			char filename[1024];

448
449
450
451
452
453
454
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

455
			number++;
Thomas White's avatar
Thomas White committed
456

457
458
			/* Write the output file */
			hdf5_write(filename, image.data,
459
			           image.width, image.height, H5T_NATIVE_FLOAT);
460
461

		}
462

463
464
465
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
466

467
skip:
Thomas White's avatar
Thomas White committed
468
		ndone++;
469

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

472
	} while ( !done );
473

474
475
476
477
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
478
479
	free(image.det->panels);
	free(image.det);
480
	free(powder);
481
	cell_free(cell);
482
	free(intensities);
Thomas White's avatar
Thomas White committed
483
484
	free(outfile);
	free(filename);
485

486
	return 0;
487
}