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


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


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

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


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

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

217
	/* Short options */
218
	while ((c = getopt_long(argc, argv, "hrn:i:g:p:o:", longopts, NULL)) != -1) {
219

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

226
227
228
229
230
231
232
233
234
235
		case 'r' : {
			config_randomquat = 1;
			break;
		}

		case 'n' : {
			n_images = atoi(optarg);
			break;
		}

236
237
238
239
240
		case 'i' : {
			intfile = strdup(optarg);
			break;
		}

241
242
243
244
245
		case 'g' : {
			grad_str = strdup(optarg);
			break;
		}

246
247
248
249
250
		case 'p' : {
			filename = strdup(optarg);
			break;
		}

251
252
253
254
255
		case 'o' : {
			outfile = strdup(optarg);
			break;
		}

256
257
258
259
260
		case 'w' : {
			powder_fn = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
261
262
263
		case 0 : {
			break;
		}
264

Thomas White's avatar
Thomas White committed
265
266
267
		default : {
			return 1;
		}
268
269
270
271
		}

	}

272
273
274
275
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

276
277
278
279
280
281
282
283
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

Thomas White's avatar
Thomas White committed
284
	if ( config_simdetails ) {
285
		show_details();
Thomas White's avatar
Thomas White committed
286
287
288
		return 0;
	}

289
290
291
292
293
294
	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;
	}

295
296
297
298
299
300
	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;
301
		counts = NULL;
302
		phases = NULL;
303
	} else {
304
		counts = new_list_count();
305
306
		phases = new_list_phase();
		intensities = read_reflections(intfile, counts, phases);
307
308
309
		free(intfile);
	}

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
	/* Define image parameters */
Thomas White's avatar
Thomas White committed
333
334
	image.width = 1024;
	image.height = 1024;
335
	image.lambda = ph_en_to_lambda(eV_to_J(PHOTON_ENERGY)); /* Wavelength */
336
	cell = load_cell_from_pdb(filename);
Thomas White's avatar
Thomas White committed
337
338
339
	if ( cell == NULL ) {
		exit(1);
	}
340
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
341
342
	image.features = NULL;
	image.flags = NULL;
343
344
	image.f0 = 1.0;
	image.f0_available = 1;
Thomas White's avatar
Thomas White committed
345

Thomas White's avatar
Thomas White committed
346
	#include "geometry-lcls.tmp"
347

348
349
	powder = calloc(image.width*image.height, sizeof(*powder));

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

353
354
	do {

Thomas White's avatar
Thomas White committed
355
		int na, nb, nc;
356
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
357

358
359
360
		//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
361
362
363
		na = 24;
		nb = 24;
		nc = 40;
Thomas White's avatar
Thomas White committed
364

365
366
367
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
			image.orientation = random_quaternion();
368
		} else {
369
			image.orientation = read_quaternion();
370
371
		}

Thomas White's avatar
Thomas White committed
372
373
374
375
		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
376
		if ( !quaternion_valid(image.orientation) ) {
377
			ERROR("Orientation modulus is not zero!\n");
378
379
380
381
382
383
			return 1;
		}

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

385
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
386
387
388
		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);

389
		if ( config_gpu ) {
390
391
			if ( gctx == NULL ) {
				gctx = setup_gpu(config_nosfac, &image,
392
				                 intensities, counts);
393
			}
394
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
395
		} else {
396
			get_diffraction(&image, na, nb, nc, intensities, counts,
397
			                phases, cell, !config_nowater, grad);
398
		}
399
		if ( image.data == NULL ) {
400
401
402
403
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

404
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
405

Thomas White's avatar
Thomas White committed
406
		if ( config_nearbragg ) {
Thomas White's avatar
Thomas White committed
407
			find_projected_peaks(&image, cell);
408
			output_intensities(&image, cell, NULL, 1);
Thomas White's avatar
Thomas White committed
409
410
		}

411
		if ( powder_fn != NULL ) {
412
413
414
415
416
417
418

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
419
				powder[x+w*y] += (double)image.data[x+w*y];
420
421
422
423
			}
			}

			if ( !(ndone % 10) ) {
424
				hdf5_write(powder_fn, powder,
425
				           image.width, image.height,
426
				           H5T_NATIVE_DOUBLE);
427
428
429
			}
		}

430
431
432
433
		if ( !config_noimages ) {

			char filename[1024];

434
435
436
437
438
439
440
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

441
			number++;
Thomas White's avatar
Thomas White committed
442

443
444
			/* Write the output file */
			hdf5_write(filename, image.data,
445
			           image.width, image.height, H5T_NATIVE_FLOAT);
446
447

		}
448

449
450
451
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
452

453
skip:
Thomas White's avatar
Thomas White committed
454
		ndone++;
455

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

458
	} while ( !done );
459

460
461
462
463
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

464
465
	free(image.det.panels);
	free(powder);
466
467
	free(cell);
	free(intensities);
468

469
	return 0;
470
}