pattern_sim.c 11.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 "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
65
" -i, --intensities=<file>  Specify file containing reflection intensities\n"
"                            to use.\n"
66
67
68
69
70
71
72
73
" -g, --gradients=<method>  Use <method> for the calculation of shape\n"
"                            transform intensities.  Choose from:\n"
"                             mosaic      : Take the intensity of the nearest\n"
"                                           Bragg position.\n"
"                             interpolate : Interpolate trilinearly between\n"
"                                           six adjacent Bragg intensities.\n"
"                             phased      : As 'interpolate', but take phase\n"
"                                           values into account.\n"
Thomas White's avatar
Thomas White committed
74
75
76
77
78
"\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"
79
"     --no-water            Do not simulate water background.\n"
Thomas White's avatar
Thomas White committed
80
"     --no-noise            Do not calculate Poisson noise.\n"
81
);
Thomas White's avatar
Thomas White committed
82
83
84
}


85
static void show_details()
Thomas White's avatar
Thomas White committed
86
{
Thomas White's avatar
Thomas White committed
87
88
89
90
91
92
93
94
95
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"
"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"
101
"This square modulus of this value is taken, and multiplied by the Bragg\n"
102
103
104
105
"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
106
"\n"
107
108
109
"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
110
111
112
113
114
115
116
117
"\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"
118
"where |F(q)|^2 is the value calculated as described above.\n"
Thomas White's avatar
Thomas White committed
119
120
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
121
122
123
"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"
124
);
125
126
127
128
129
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
}


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 {
155
			ERROR("Invalid rotation '%s'\n", line);
156
157
158
		}

	} while ( 1 );
159
160
161
162
163
}


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

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

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

Thomas White's avatar
Thomas White committed
212
213
214
215
216
		switch (c) {
		case 'h' : {
			show_help(argv[0]);
			return 0;
		}
217

218
219
220
221
222
223
224
225
226
227
		case 'r' : {
			config_randomquat = 1;
			break;
		}

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

228
229
230
231
232
		case 'i' : {
			intfile = strdup(optarg);
			break;
		}

233
234
235
236
237
		case 'g' : {
			grad_str = strdup(optarg);
			break;
		}

238
239
240
241
242
		case 'p' : {
			filename = strdup(optarg);
			break;
		}

243
244
245
246
247
		case 'o' : {
			outfile = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
248
249
250
		case 0 : {
			break;
		}
251

Thomas White's avatar
Thomas White committed
252
253
254
		default : {
			return 1;
		}
255
256
257
258
		}

	}

259
260
261
262
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

263
264
265
266
267
268
269
270
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

Thomas White's avatar
Thomas White committed
271
	if ( config_simdetails ) {
272
		show_details();
Thomas White's avatar
Thomas White committed
273
274
275
		return 0;
	}

276
277
278
279
280
281
	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;
	}

282
283
284
285
286
287
	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;
288
		counts = NULL;
289
	} else {
290
291
		counts = new_list_count();
		intensities = read_reflections(intfile, counts);
292
293
294
		free(intfile);
	}

295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
	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
317
	/* Define image parameters */
Thomas White's avatar
Thomas White committed
318
319
	image.width = 1024;
	image.height = 1024;
320
	image.lambda = ph_en_to_lambda(eV_to_J(PHOTON_ENERGY)); /* Wavelength */
321
	cell = load_cell_from_pdb(filename);
322
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
323
324
	image.features = NULL;
	image.flags = NULL;
Thomas White's avatar
Thomas White committed
325

Thomas White's avatar
Thomas White committed
326
	#include "geometry-lcls.tmp"
327

328
329
	powder = calloc(image.width*image.height, sizeof(*powder));

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

333
334
	do {

Thomas White's avatar
Thomas White committed
335
		int na, nb, nc;
336
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
337

338
339
340
		//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
341
342
343
		na = 24;
		nb = 24;
		nc = 40;
Thomas White's avatar
Thomas White committed
344

345
346
347
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
			image.orientation = random_quaternion();
348
		} else {
349
			image.orientation = read_quaternion();
350
351
		}

Thomas White's avatar
Thomas White committed
352
353
354
355
		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
356
		if ( !quaternion_valid(image.orientation) ) {
357
			ERROR("Orientation modulus is not zero!\n");
358
359
360
361
362
363
			return 1;
		}

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

365
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
366
367
368
		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);

369
		if ( config_gpu ) {
370
371
			if ( gctx == NULL ) {
				gctx = setup_gpu(config_nosfac, &image,
372
				                 intensities, counts);
373
			}
374
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
375
		} else {
376
			get_diffraction(&image, na, nb, nc, intensities, counts,
377
			                cell, !config_nowater, grad);
378
		}
379
		if ( image.data == NULL ) {
380
381
382
383
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

384
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
385

Thomas White's avatar
Thomas White committed
386
		if ( config_nearbragg ) {
Thomas White's avatar
Thomas White committed
387
			find_projected_peaks(&image, cell);
388
			output_intensities(&image, cell, NULL, 1);
Thomas White's avatar
Thomas White committed
389
390
		}

391
392
393
394
395
396
397
398
		if ( config_powder ) {

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
399
				powder[x+w*y] += (double)image.data[x+w*y];
400
401
402
403
			}
			}

			if ( !(ndone % 10) ) {
404
				hdf5_write("results/integr.h5", powder,
405
				           image.width, image.height,
406
				           H5T_NATIVE_DOUBLE);
407
408
409
			}
		}

410
411
412
413
		if ( !config_noimages ) {

			char filename[1024];

414
415
416
417
418
419
420
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

421
			number++;
Thomas White's avatar
Thomas White committed
422

423
424
			/* Write the output file */
			hdf5_write(filename, image.data,
425
			           image.width, image.height, H5T_NATIVE_FLOAT);
426
427

		}
428

429
430
431
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
432

433
skip:
Thomas White's avatar
Thomas White committed
434
		ndone++;
435

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

438
	} while ( !done );
439

440
441
442
443
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

444
445
	free(image.det.panels);
	free(powder);
446
447
	free(cell);
	free(intensities);
448

449
	return 0;
450
}