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"
Thomas White's avatar
Thomas White committed
62
63
"     --powder              Write a summed pattern of all images simulated by\n"
"                            this invocation to results/integr.h5.\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
98
99
100
101
102
103
104
105
106
	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"
"using the closed-form solution for a truncated lattice:\n"
"\n"
"F(q) =  sin(pi*na*q.a)/sin(pi*q.a)\n"
"      * sin(pi*nb*q.b)/sin(pi*q.b)\n"
"      * sin(pi*nc*q.c)/sin(pi*q.c)\n"
"\n"
"na = number of unit cells in 'a' direction (likewise nb, nc)\n"
" q = reciprocal vector (1/d convention, not 2pi/d)\n"
"\n"
107
"This square modulus of this value is taken, and multiplied by the Bragg\n"
108
109
110
111
"intensitiy at the neares Bragg position.  That means that the gradient of\n"
"the underlying molecule transform is not included, for speed of calculation.\n"
"The Bragg intensities are taken from the file specified on the command line\n"
"with the --intensities option.\n"
Thomas White's avatar
Thomas White committed
112
"\n"
113
114
115
"Intensity from water is then added according to the first term of equation\n"
"5 from Phys Chem Chem Phys 2003 (5) 1981--1991.  This simulates the\n"
"coherent, elastic part of the diffuse scattering from the water jet only.\n"
Thomas White's avatar
Thomas White committed
116
117
118
119
120
121
122
123
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
"I(q) = I0 * r^2 * |F(q)|^2 * S\n"
"\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
"where |F(q)|^2 is the value calculated as described above.\n"
Thomas White's avatar
Thomas White committed
125
126
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
127
128
129
"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"
130
);
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
}


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

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


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

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

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

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

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

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

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

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

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

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

Thomas White's avatar
Thomas White committed
255
256
257
		case 0 : {
			break;
		}
258

Thomas White's avatar
Thomas White committed
259
260
261
		default : {
			return 1;
		}
262
263
264
265
		}

	}

266
267
268
269
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

270
271
272
273
274
275
276
277
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

Thomas White's avatar
Thomas White committed
278
	if ( config_simdetails ) {
279
		show_details();
Thomas White's avatar
Thomas White committed
280
281
282
		return 0;
	}

283
284
285
286
287
288
	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;
	}

289
290
291
292
293
294
	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;
295
		counts = NULL;
296
		phases = NULL;
297
	} else {
298
		counts = new_list_count();
299
300
		phases = new_list_phase();
		intensities = read_reflections(intfile, counts, phases);
301
302
303
		free(intfile);
	}

304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
	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
326
	/* Define image parameters */
Thomas White's avatar
Thomas White committed
327
328
	image.width = 1024;
	image.height = 1024;
329
	image.lambda = ph_en_to_lambda(eV_to_J(PHOTON_ENERGY)); /* Wavelength */
330
	cell = load_cell_from_pdb(filename);
Thomas White's avatar
Thomas White committed
331
332
333
	if ( cell == NULL ) {
		exit(1);
	}
334
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
335
336
	image.features = NULL;
	image.flags = NULL;
337
338
	image.f0 = 1.0;
	image.f0_available = 1;
Thomas White's avatar
Thomas White committed
339

Thomas White's avatar
Thomas White committed
340
	#include "geometry-lcls.tmp"
341

342
343
	powder = calloc(image.width*image.height, sizeof(*powder));

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

347
348
	do {

Thomas White's avatar
Thomas White committed
349
		int na, nb, nc;
350
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
351

352
353
354
		//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
355
356
357
		na = 24;
		nb = 24;
		nc = 40;
Thomas White's avatar
Thomas White committed
358

359
360
361
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
			image.orientation = random_quaternion();
362
		} else {
363
			image.orientation = read_quaternion();
364
365
		}

Thomas White's avatar
Thomas White committed
366
367
368
369
		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
370
		if ( !quaternion_valid(image.orientation) ) {
371
			ERROR("Orientation modulus is not zero!\n");
372
373
374
375
376
377
			return 1;
		}

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

379
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
380
381
382
		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);

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

398
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
399

Thomas White's avatar
Thomas White committed
400
		if ( config_nearbragg ) {
Thomas White's avatar
Thomas White committed
401
			find_projected_peaks(&image, cell);
402
			output_intensities(&image, cell, NULL, 1);
Thomas White's avatar
Thomas White committed
403
404
		}

405
406
407
408
409
410
411
412
		if ( config_powder ) {

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
413
				powder[x+w*y] += (double)image.data[x+w*y];
414
415
416
417
			}
			}

			if ( !(ndone % 10) ) {
418
				hdf5_write("results/integr.h5", powder,
419
				           image.width, image.height,
420
				           H5T_NATIVE_DOUBLE);
421
422
423
			}
		}

424
425
426
427
		if ( !config_noimages ) {

			char filename[1024];

428
429
430
431
432
433
434
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

435
			number++;
Thomas White's avatar
Thomas White committed
436

437
438
			/* Write the output file */
			hdf5_write(filename, image.data,
439
			           image.width, image.height, H5T_NATIVE_FLOAT);
440
441

		}
442

443
444
445
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
446

447
skip:
Thomas White's avatar
Thomas White committed
448
		ndone++;
449

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

452
	} while ( !done );
453

454
455
456
457
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

458
459
	free(image.det.panels);
	free(powder);
460
461
	free(cell);
	free(intensities);
462

463
	return 0;
464
}