pattern_sim.c 10.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
35


Thomas White's avatar
Thomas White committed
36
static void show_help(const char *s)
37
{
38
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
39
	printf(
Thomas White's avatar
Thomas White committed
40
"Simulate diffraction patterns from small crystals probed with femtosecond\n"
Thomas White's avatar
Thomas White committed
41
42
"pulses of X-rays from a free electron laser.\n"
"\n"
43
" -h, --help                Display this help message.\n"
44
45
46
" -p, --pdb=<file>          PDB file from which to get the unit cell.\n"
"                            (The actual Bragg intensities come from the\n"
"                            intensities file)\n"
47
"     --simulation-details  Show technical details of the simulation.\n"
48
"     --gpu                 Use the GPU to speed up the calculation.\n"
49
50
51
52
"\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"
Thomas White's avatar
Thomas White committed
53
" -r, --random-orientation  Use a randomly generated orientation\n"
54
"                            (a new orientation will be used for each image).\n"
Thomas White's avatar
Thomas White committed
55
56
"     --powder              Write a summed pattern of all images simulated by\n"
"                            this invocation to results/integr.h5.\n"
57
58
" -i, --intensities=<file>  Specify file containing reflection intensities\n"
"                            to use.\n"
59
60
61
62
63
64
65
66
" -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
67
68
69
70
71
"\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"
72
"     --no-water            Do not simulate water background.\n"
Thomas White's avatar
Thomas White committed
73
"     --no-noise            Do not calculate Poisson noise.\n"
74
);
Thomas White's avatar
Thomas White committed
75
76
77
}


78
static void show_details()
Thomas White's avatar
Thomas White committed
79
{
Thomas White's avatar
Thomas White committed
80
81
82
83
84
85
86
87
88
89
90
91
92
93
	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"
94
"This square modulus of this value is taken, and multiplied by the Bragg\n"
95
96
97
98
"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
99
"\n"
100
101
102
"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
103
104
105
106
107
108
109
110
"\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"
111
"where |F(q)|^2 is the value calculated as described above.\n"
Thomas White's avatar
Thomas White committed
112
113
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
114
115
116
"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"
117
);
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
}


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 {
148
			ERROR("Invalid rotation '%s'\n", line);
149
150
151
		}

	} while ( 1 );
152
153
154
155
156
}


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

182
	/* Long options */
Thomas White's avatar
Thomas White committed
183
	const struct option longopts[] = {
184
185
		{"help",               0, NULL,               'h'},
		{"simulation-details", 0, &config_simdetails,  1},
186
		{"gpu",                0, &config_gpu,         1},
187
188
189
		{"near-bragg",         0, &config_nearbragg,   1},
		{"random-orientation", 0, NULL,               'r'},
		{"number",             1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
190
		{"no-images",          0, &config_noimages,    1},
191
		{"no-water",           0, &config_nowater,     1},
Thomas White's avatar
Thomas White committed
192
		{"no-noise",           0, &config_nonoise,     1},
193
		{"intensities",        1, NULL,               'i'},
194
		{"powder",             0, &config_powder,      1},
195
		{"gradients",          1, NULL,               'g'},
196
		{"pdb",                1, NULL,               'p'},
197
		{0, 0, NULL, 0}
Thomas White's avatar
Thomas White committed
198
	};
199

200
	/* Short options */
201
	while ((c = getopt_long(argc, argv, "hrn:i:g:p:", longopts, NULL)) != -1) {
202

Thomas White's avatar
Thomas White committed
203
204
205
206
207
		switch (c) {
		case 'h' : {
			show_help(argv[0]);
			return 0;
		}
208

209
210
211
212
213
214
215
216
217
218
		case 'r' : {
			config_randomquat = 1;
			break;
		}

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

219
220
221
222
223
		case 'i' : {
			intfile = strdup(optarg);
			break;
		}

224
225
226
227
228
		case 'g' : {
			grad_str = strdup(optarg);
			break;
		}

229
230
231
232
233
		case 'p' : {
			filename = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
234
235
236
		case 0 : {
			break;
		}
237

Thomas White's avatar
Thomas White committed
238
239
240
		default : {
			return 1;
		}
241
242
243
244
		}

	}

245
246
247
248
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

Thomas White's avatar
Thomas White committed
249
	if ( config_simdetails ) {
250
		show_details();
Thomas White's avatar
Thomas White committed
251
252
253
		return 0;
	}

254
255
256
257
258
259
	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;
	}

260
261
262
263
264
265
	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;
266
		counts = NULL;
267
	} else {
268
269
		counts = new_list_count();
		intensities = read_reflections(intfile, counts);
270
271
272
		free(intfile);
	}

273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
	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
295
	/* Define image parameters */
Thomas White's avatar
Thomas White committed
296
297
	image.width = 1024;
	image.height = 1024;
Thomas White's avatar
Thomas White committed
298
	image.lambda = ph_en_to_lambda(eV_to_J(1790.0));  /* Wavelength */
299
	cell = load_cell_from_pdb(filename);
300
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
301

Thomas White's avatar
Thomas White committed
302
	#include "geometry-lcls.tmp"
303

304
305
	powder = calloc(image.width*image.height, sizeof(*powder));

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

309
310
	do {

Thomas White's avatar
Thomas White committed
311
		int na, nb, nc;
312
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
313

314
315
316
		//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
317
318
319
		na = 24;
		nb = 24;
		nc = 40;
Thomas White's avatar
Thomas White committed
320

321
322
323
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
			image.orientation = random_quaternion();
324
		} else {
325
			image.orientation = read_quaternion();
326
327
		}

Thomas White's avatar
Thomas White committed
328
329
330
331
		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
332
		if ( !quaternion_valid(image.orientation) ) {
333
			ERROR("Orientation modulus is not zero!\n");
334
335
336
337
338
339
			return 1;
		}

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

341
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
342
343
344
		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);

345
		if ( config_gpu ) {
346
347
			if ( gctx == NULL ) {
				gctx = setup_gpu(config_nosfac, &image,
348
				                 intensities, counts);
349
			}
350
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
351
		} else {
352
			get_diffraction(&image, na, nb, nc, intensities, counts,
353
			                cell, !config_nowater, grad);
354
		}
355
		if ( image.data == NULL ) {
356
357
358
359
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

360
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
361

Thomas White's avatar
Thomas White committed
362
		if ( config_nearbragg ) {
363
			output_intensities(&image, cell, NULL, 1);
Thomas White's avatar
Thomas White committed
364
365
		}

366
367
368
369
370
371
372
373
		if ( config_powder ) {

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
374
				powder[x+w*y] += (double)image.data[x+w*y];
375
376
377
378
			}
			}

			if ( !(ndone % 10) ) {
379
				hdf5_write("results/integr.h5", powder,
380
				           image.width, image.height,
381
				           H5T_NATIVE_DOUBLE);
382
383
384
			}
		}

385
386
387
388
389
390
		if ( !config_noimages ) {

			char filename[1024];

			snprintf(filename, 1023, "results/sim-%i.h5", number);
			number++;
Thomas White's avatar
Thomas White committed
391

392
393
			/* Write the output file */
			hdf5_write(filename, image.data,
394
			           image.width, image.height, H5T_NATIVE_FLOAT);
395
396

		}
397

398
399
400
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
401

402
skip:
Thomas White's avatar
Thomas White committed
403
		ndone++;
404

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

407
	} while ( !done );
408

409
410
411
412
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

413
414
	free(image.det.panels);
	free(powder);
415
416
	free(cell);
	free(intensities);
417

418
	return 0;
419
}