get_hkl.c 10.2 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
static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
111
112
                                      const char *holo, const char *mero,
                                      double *esds)
113
114
115
116
117
118
{
	int i;
	ReflItemList *new;

	new = new_items();

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

124
125
	for ( i=0; i<num_items(items); i++ ) {

126
		double total, sigma;
127
128
129
130
131
132
133
		struct refl_item *it;
		signed int h, k, l;
		int n, j;
		int skip;

		it = get_item(items, i);

134
		/* There is a many-to-one correspondence between reflections
Thomas White's avatar
Thomas White committed
135
136
137
		 * in the merohedral and holohedral groups.  Do the calculation
		 * only once for each reflection in the holohedral group, which
		 * contains fewer reflections.
138
		 */
139
140
		get_asymm(it->h, it->k, it->l, &h, &k, &l, holo);
		if ( find_item(new, h, k, l) ) continue;
141
142
143

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

144
145
		total = 0.0;
		sigma = 0.0;
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
		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;
			}

171
172
			total += lookup_intensity(ref, hu, ku, lu);
			sigma += pow(lookup_intensity(esds, hu, ku, lu), 2.0);
173
174
175
176
177

		}

		if ( !skip ) {

178
179
			set_intensity(ref, h, k, l, total);
			set_intensity(esds, h, k, l, sqrt(sigma));
180
181
182
183
184
185
186
187
188
189
			add_item(new, h, k, l);

		}

	}

	return new;
}


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

	new = new_items();

198
199
200
201
202
	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
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
247
248
	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;
}


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

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

Thomas White's avatar
Thomas White committed
325
		case 0 :
326
327
			break;

Thomas White's avatar
Thomas White committed
328
		default :
329
330
331
332
333
			return 1;
		}

	}

334
335
336
337
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

Thomas White's avatar
Thomas White committed
338
339
340
341
342
343
	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);
	}

344
	mol = load_molecule(filename);
345
	cell = load_cell_from_pdb(filename);
346
347
348
349
350
	if ( !config_nophase ) {
		phases = new_list_phase();
	} else {
		phases = NULL;
	}
351
	esds = new_list_intensity();
352
	if ( input == NULL ) {
Thomas White's avatar
Thomas White committed
353
		input_items = new_items();
354
		ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9),
Thomas White's avatar
Thomas White committed
355
		                            phases, input_items);
356
	} else {
Thomas White's avatar
Thomas White committed
357
		ideal_ref = new_list_intensity();
358
		input_items = read_reflections(input, ideal_ref, phases,
359
		                               NULL, esds);
360
		free(input);
Thomas White's avatar
Thomas White committed
361
362
363
364
365
		if ( check_symmetry(input_items, mero) ) {
			ERROR("The input reflection list does not appear to"
			      " have symmetry %s\n", mero);
			return 1;
		}
366
	}
367

368
369
	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
370

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

Thomas White's avatar
Thomas White committed
380
381
382
383
384
385
386
387
	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;
	}

388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
	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
407
408
409
410
	if ( template ) {
		/* Write out only reflections which are in the template
		 * (and which we have in the input) */
		ReflItemList *template_items;
411
412
		template_items = read_reflections(template,
		                                  NULL, NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
413
414
415
416
417
		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
418
		/* (quick way of copying a list) */
Thomas White's avatar
Thomas White committed
419
420
421
		union_items(write_items, input_items);
	}

422
	write_reflections(output, write_items, ideal_ref, esds, phases,
423
	                  NULL, cell);
Thomas White's avatar
Thomas White committed
424
425
426

	delete_items(input_items);
	delete_items(write_items);
427
428
429

	return 0;
}