get_hkl.c 10.4 KB
Newer Older
1
2
3
4
5
/*
 * get_hkl.c
 *
 * Small program to write out a list of h,k,l,I values given a structure
 *
Thomas White's avatar
Thomas White committed
6
 * (c) 2006-2010 Thomas White <taw@physics.org>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 *
 * Part of CrystFEL - crystallography with a FEL
 *
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <getopt.h>

#include "utils.h"
#include "sfac.h"
#include "reflections.h"
27
#include "symmetry.h"
28
#include "beam-parameters.h"
29
30
31
32
33
34


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
35
"Create reflections lists.\n"
36
"\n"
Thomas White's avatar
Thomas White committed
37
38
"  -h, --help                 Display this help message.\n"
"\n"
Thomas White's avatar
Thomas White committed
39
40
"  -t, --template=<filename>  Only include reflections mentioned in file.\n"
"      --poisson              Simulate Poisson samples.\n"
41
"      --noise                Add 10%% random noise.\n"
42
43
44
"  -y, --symmetry=<sym>       The symmetry of the input file (-i).\n"
"  -w, --twin=<sym>           Generate twinned data according to the given\n"
"                              point group.\n"
Thomas White's avatar
Thomas White committed
45
"  -e, --expand=<sym>         Expand reflections to this point group.\n"
46
47
"  -o, --output=<filename>    Output filename (default: stdout).\n"
"  -i, --intensities=<file>   Read intensities from file instead of\n"
48
49
50
"                              calculating them from scratch.  You might use\n"
"                              this if you need to apply noise or twinning.\n"
"  -p, --pdb=<file>           PDB file from which to get the structure.\n"
51
52
53
"      --no-phases            Do not try to use phases in the input file.\n"
"      --multiplicity         Multiply intensities by the number of\n"
"                              equivalent reflections.\n"
54
"  -b, --beam=<file>         Get beam parameters from file (used for sigmas).\n"
Thomas White's avatar
Thomas White committed
55
);
Thomas White's avatar
Thomas White committed
56
57
58
59
}


/* Apply Poisson noise to all reflections */
60
static void poisson_reflections(double *ref, ReflItemList *items)
Thomas White's avatar
Thomas White committed
61
{
62
63
	int i;
	const int n = num_items(items);
Thomas White's avatar
Thomas White committed
64

65
	for ( i=0; i<n; i++ ) {
Thomas White's avatar
Thomas White committed
66

67
		struct refl_item *it;
Thomas White's avatar
Thomas White committed
68
		double val;
Thomas White's avatar
Thomas White committed
69
		int c;
Thomas White's avatar
Thomas White committed
70

71
72
73
		it = get_item(items, i);

		val = lookup_intensity(ref, it->h, it->k, it->l);
Thomas White's avatar
Thomas White committed
74
		c = poisson_noise(val);
75
76
77
		set_intensity(ref, it->h, it->k, it->l, c);

		progress_bar(i, n-1, "Simulating noise");
Thomas White's avatar
Thomas White committed
78
79

	}
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
}


/* Apply 10% uniform noise to all reflections */
static void noise_reflections(double *ref, ReflItemList *items)
{
	int i;
	const int n = num_items(items);

	for ( i=0; i<n; i++ ) {

		struct refl_item *it;
		double val;
		double r;

		it = get_item(items, i);

		val = lookup_intensity(ref, it->h, it->k, it->l);

		r = (double)random()/RAND_MAX;
		val += 0.1 * val * r;

		set_intensity(ref, it->h, it->k, it->l, val);

		progress_bar(i, n-1, "Simulating noise");

Thomas White's avatar
Thomas White committed
106
	}
107
108
109
}


110
111
112
113
114
115
116
117
static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
                                      const char *holo, const char *mero)
{
	int i;
	ReflItemList *new;

	new = new_items();

118
119
120
121
122
	if ( num_general_equivs(holo) < num_general_equivs(mero) ) {
		ERROR("%s is not a subgroup of %s!\n", mero, holo);
		return NULL;
	}

123
124
125
126
127
128
129
130
131
132
	for ( i=0; i<num_items(items); i++ ) {

		double mean;
		struct refl_item *it;
		signed int h, k, l;
		int n, j;
		int skip;

		it = get_item(items, i);

133
		/* There is a many-to-one correspondence between reflections
Thomas White's avatar
Thomas White committed
134
135
136
		 * in the merohedral and holohedral groups.  Do the calculation
		 * only once for each reflection in the holohedral group, which
		 * contains fewer reflections.
137
		 */
138
139
		get_asymm(it->h, it->k, it->l, &h, &k, &l, holo);
		if ( find_item(new, h, k, l) ) continue;
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187

		n = num_equivs(h, k, l, holo);

		mean = 0.0;
		skip = 0;
		for ( j=0; j<n; j++ ) {

			signed int he, ke, le;
			signed int hu, ku, lu;

			get_equiv(h, k, l, &he, &ke, &le, holo, j);

			/* Do we have this reflection?
			 * We might not have the particular (merohedral)
			 * equivalent which belongs to our definition of the
			 * asymmetric unit cell, so check them all.
			 */
			if ( !find_unique_equiv(items, he, ke, le, mero,
			                        &hu, &ku, &lu) ) {
				/* Don't have this reflection, so bail out */
				ERROR("Twinning %i %i %i requires the %i %i %i "
				      "reflection (or an equivalent in %s), "
				      "which I don't have. %i %i %i won't "
				      "appear in the output\n",
				      h, k, l, he, ke, le, mero, h, k, l);
				skip = 1;
				break;
			}

			mean += lookup_intensity(ref, hu, ku, lu);

		}

		if ( !skip ) {

			mean /= (double)n;

			set_intensity(ref, h, k, l, mean);
			add_item(new, h, k, l);

		}

	}

	return new;
}


Thomas White's avatar
Thomas White committed
188
189
190
191
192
193
194
195
static ReflItemList *expand_reflections(double *ref, ReflItemList *items,
                                        const char *target, const char *initial)
{
	int i;
	ReflItemList *new;

	new = new_items();

196
197
198
199
200
	if ( num_general_equivs(target) > num_general_equivs(initial) ) {
		ERROR("%s is not a subgroup of %s!\n", initial, target);
		return NULL;
	}

Thomas White's avatar
Thomas White committed
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
	for ( i=0; i<num_items(items); i++ ) {

		struct refl_item *it;
		signed int h, k, l;
		signed int hd, kd, ld;
		int n, j;
		double intensity;

		it = get_item(items, i);
		h = it->h;  k = it->k;  l = it->l;

		/* Actually we don't really care what the equivalent is,
		 * we just want to be sure that there is nly be one version of
		 * this reflection. */
		find_unique_equiv(items, h, k, l, initial, &hd, &kd, &ld);

		/* Now find out how many reflections need to be filled in */
		n = num_equivs(h, k, l, initial);
		intensity = lookup_intensity(ref, h, k, l);

		for ( j=0; j<n; j++ ) {

			signed int he, ke, le;

			/* Get the equivalent */
			get_equiv(h, k, l, &he, &ke, &le, initial, j);

			/* Put it into the asymmetric unit for the target */
			get_asymm(he, ke, le, &he, &ke, &le, target);

			/* Make sure the intensity is in the right place */
			set_intensity(ref, he, ke, le, intensity);

			/* Add the reflection, but only once */
			if ( !find_item(new, he, ke, le) ) {
				add_item(new, he, ke, le);
			}

		}

	}

	return new;
}


247
248
249
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
250
	double *ideal_ref;
251
	double *phases;
252
	struct molecule *mol;
Thomas White's avatar
Thomas White committed
253
	char *template = NULL;
254
255
	int config_noise = 0;
	int config_poisson = 0;
256
257
	int config_nophase = 0;
	int config_multi = 0;
258
259
	char *holo = NULL;
	char *mero = NULL;
Thomas White's avatar
Thomas White committed
260
	char *expand = NULL;
Thomas White's avatar
Thomas White committed
261
	char *output = NULL;
262
	char *input = NULL;
263
	char *filename = NULL;
Thomas White's avatar
Thomas White committed
264
265
	ReflItemList *input_items;
	ReflItemList *write_items;
266
	UnitCell *cell = NULL;
267
268
	double adu_per_photon;
	struct beam_params *beam = NULL;
269
270
271
272

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
273
		{"template",           1, NULL,               't'},
274
275
		{"poisson",            0, &config_poisson,     1},
		{"noise",              0, &config_noise,       1},
Thomas White's avatar
Thomas White committed
276
		{"output",             1, NULL,               'o'},
277
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
278
279
		{"twin",               1, NULL,               'w'},
		{"expand",             1, NULL,               'e'},
280
		{"intensities",        1, NULL,               'i'},
281
		{"pdb",                1, NULL,               'p'},
282
283
		{"no-phases",          0, &config_nophase,     1},
		{"multiplicity",       0, &config_multi,       1},
284
		{"beam",               1, NULL,               'b'},
285
286
287
288
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
289
	while ((c = getopt_long(argc, argv, "ht:o:i:p:w:y:e:",
Thomas White's avatar
Thomas White committed
290
	                        longopts, NULL)) != -1) {
291
292

		switch (c) {
Thomas White's avatar
Thomas White committed
293
		case 'h' :
294
295
296
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
297
		case 't' :
Thomas White's avatar
Thomas White committed
298
299
300
			template = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
301
		case 'o' :
Thomas White's avatar
Thomas White committed
302
303
304
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
305
		case 'i' :
306
307
308
			input = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
309
		case 'p' :
310
311
312
			filename = strdup(optarg);
			break;

313
314
315
316
317
318
319
320
		case 'y' :
			mero = strdup(optarg);
			break;

		case 'w' :
			holo = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
321
322
323
324
		case 'e' :
			expand = strdup(optarg);
			break;

325
326
327
328
329
330
331
332
333
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
334
		case 0 :
335
336
			break;

Thomas White's avatar
Thomas White committed
337
		default :
338
339
340
341
342
			return 1;
		}

	}

343
344
345
346
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

Thomas White's avatar
Thomas White committed
347
348
349
350
351
352
	if ( (holo != NULL) && (expand != NULL) ) {
		ERROR("You cannot 'twin' and 'expand' at the same time.\n");
		ERROR("Decide which one you want to do first.\n");
		exit(1);
	}

353
	mol = load_molecule(filename);
354
	cell = load_cell_from_pdb(filename);
355
356
357
358
359
	if ( !config_nophase ) {
		phases = new_list_phase();
	} else {
		phases = NULL;
	}
360
	if ( input == NULL ) {
Thomas White's avatar
Thomas White committed
361
		input_items = new_items();
362
		ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9),
Thomas White's avatar
Thomas White committed
363
		                            phases, input_items);
364
	} else {
Thomas White's avatar
Thomas White committed
365
		ideal_ref = new_list_intensity();
366
367
		input_items = read_reflections(input, ideal_ref, phases,
		                               NULL, NULL);
368
369
		free(input);
	}
370

371
372
	if ( config_poisson ) poisson_reflections(ideal_ref, input_items);
	if ( config_noise ) noise_reflections(ideal_ref, input_items);
Thomas White's avatar
Thomas White committed
373

374
375
376
377
378
379
	if ( holo != NULL ) {
		ReflItemList *new;
		STATUS("Twinning from %s into %s\n", mero, holo);
		new = twin_reflections(ideal_ref, input_items, holo, mero);
		delete_items(input_items);
		input_items = new;
Thomas White's avatar
Thomas White committed
380
381
	}

Thomas White's avatar
Thomas White committed
382
383
384
385
386
387
388
389
	if ( expand != NULL ) {
		ReflItemList *new;
		STATUS("Expanding from %s into %s\n", mero, expand);
		new = expand_reflections(ideal_ref, input_items, expand, mero);
		delete_items(input_items);
		input_items = new;
	}

390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
	if ( config_multi ) {

		int i;

		for ( i=0; i<num_items(input_items); i++ ) {

			struct refl_item *it;
			double inty;

			it = get_item(input_items, i);
			inty = lookup_intensity(ideal_ref, it->h, it->k, it->l);
			inty *= num_equivs(it->h, it->k, it->l, mero);
			set_intensity(ideal_ref, it->h, it->k, it->l, inty);
			STATUS("%i %i %i %i\n", it->h, it->k, it->l,
			       num_equivs(it->h, it->k, it->l, mero));

		}
	}

Thomas White's avatar
Thomas White committed
409
410
411
412
	if ( template ) {
		/* Write out only reflections which are in the template
		 * (and which we have in the input) */
		ReflItemList *template_items;
413
414
		template_items = read_reflections(template,
		                                  NULL, NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
415
416
417
418
419
		write_items = intersection_items(input_items, template_items);
		delete_items(template_items);
	} else {
		/* Write out all reflections */
		write_items = new_items();
Thomas White's avatar
Thomas White committed
420
		/* (quick way of copying a list) */
Thomas White's avatar
Thomas White committed
421
422
423
		union_items(write_items, input_items);
	}

424
425
426
427
428
429
430
431
	if ( beam == NULL ) {
		adu_per_photon = 167.0;
		STATUS("No beam parameters file provided (use -b), "
		       "so I'm assuming 167.0 ADU per photon.\n");
	} else {
		adu_per_photon = beam->adu_per_photon;
	}

432
433
	write_reflections(output, write_items, ideal_ref, NULL, phases,
	                  NULL, cell);
Thomas White's avatar
Thomas White committed
434
435
436

	delete_items(input_items);
	delete_items(write_items);
437
438
439

	return 0;
}