get_hkl.c 11.3 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
55
"  -b, --beam=<file>          Get beam parameters from file (used for sigmas).\n"
"      --max-res=<d>          Calculate structure factors out to d=<d> nm.\n"
Thomas White's avatar
Thomas White committed
56
);
Thomas White's avatar
Thomas White committed
57
58
59
60
}


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

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

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

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

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

		progress_bar(i, n-1, "Simulating noise");
Thomas White's avatar
Thomas White committed
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
106
}


/* 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
107
	}
108
109
110
}


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

	new = new_items();

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

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

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

		it = get_item(items, i);

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

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

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

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

		}

		if ( !skip ) {

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

		}

	}

	return new;
}


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

	new = new_items();

199
200
201
202
203
	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
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
249
	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;
}


250
251
252
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
253
	double *ideal_ref;
254
	double *phases;
255
	double *esds;
256
	struct molecule *mol;
Thomas White's avatar
Thomas White committed
257
	char *template = NULL;
258
259
	int config_noise = 0;
	int config_poisson = 0;
260
261
	int config_nophase = 0;
	int config_multi = 0;
262
263
	char *holo = NULL;
	char *mero = NULL;
Thomas White's avatar
Thomas White committed
264
	char *expand = NULL;
Thomas White's avatar
Thomas White committed
265
	char *output = NULL;
266
	char *input = NULL;
267
	char *filename = NULL;
Thomas White's avatar
Thomas White committed
268
269
	ReflItemList *input_items;
	ReflItemList *write_items;
270
	UnitCell *cell = NULL;
271
272
273
274
275
	char *beamfile = NULL;
	char *rval;
	struct beam_params *beam;  /* Beam parameters for SF calculation */
	int have_max_res = 0;
	double max_res = 0.0;
276
277
278
279

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
280
		{"template",           1, NULL,               't'},
281
282
		{"poisson",            0, &config_poisson,     1},
		{"noise",              0, &config_noise,       1},
Thomas White's avatar
Thomas White committed
283
		{"output",             1, NULL,               'o'},
284
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
285
286
		{"twin",               1, NULL,               'w'},
		{"expand",             1, NULL,               'e'},
287
		{"intensities",        1, NULL,               'i'},
288
		{"pdb",                1, NULL,               'p'},
289
290
		{"no-phases",          0, &config_nophase,     1},
		{"multiplicity",       0, &config_multi,       1},
291
292
		{"beam",               1, NULL,               'b'},
		{"max-res",            1, NULL,                2},
293
294
295
296
		{0, 0, NULL, 0}
	};

	/* Short options */
297
	while ((c = getopt_long(argc, argv, "ht:o:i:p:w:y:e:b:",
Thomas White's avatar
Thomas White committed
298
	                        longopts, NULL)) != -1) {
299
300

		switch (c) {
Thomas White's avatar
Thomas White committed
301
		case 'h' :
302
303
304
			show_help(argv[0]);
			return 0;

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

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

Thomas White's avatar
Thomas White committed
313
		case 'i' :
314
315
316
			input = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
317
		case 'p' :
318
319
320
			filename = strdup(optarg);
			break;

321
322
323
324
325
326
327
328
		case 'y' :
			mero = strdup(optarg);
			break;

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

Thomas White's avatar
Thomas White committed
329
330
331
332
		case 'e' :
			expand = strdup(optarg);
			break;

333
334
335
336
337
338
339
340
341
342
343
344
345
346
		case 'b' :
			beamfile = strdup(optarg);
			break;

		case 2 :
			max_res = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid maximum resolution.\n");
				return 1;
			}
			max_res = 1.0 / (max_res * 1.0e-9);
			have_max_res = 1;
			break;

Thomas White's avatar
Thomas White committed
347
		case 0 :
348
349
			break;

Thomas White's avatar
Thomas White committed
350
		default :
351
352
353
354
355
			return 1;
		}

	}

356
357
358
359
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

Thomas White's avatar
Thomas White committed
360
361
362
363
364
365
	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);
	}

366
	mol = load_molecule(filename);
367
	cell = load_cell_from_pdb(filename);
368
369
370
371
372
	if ( !config_nophase ) {
		phases = new_list_phase();
	} else {
		phases = NULL;
	}
373
	esds = new_list_sigma();
374
	if ( input == NULL ) {
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395

		if ( beamfile == NULL ) {
			ERROR("To calculate structure factors, you must"
			      " provide a beam parameters file (use -b)\n");
			return 1;
		}

		beam = get_beam_parameters(beamfile);
		if ( beam == NULL ) {
			ERROR("Failed to read beam parameters from '%s'\n", beamfile);
			return 1;
		}
		free(beamfile);

		if ( !have_max_res ) {
			STATUS("You didn't specify the maximum resolution to"
			       " calculate structure factors.  I'll go to"
			       " d = 0.5 nm.\n");
			max_res = 1.0/0.5e-9;
		}

Thomas White's avatar
Thomas White committed
396
		input_items = new_items();
397
398
399
		ideal_ref = get_reflections(mol, eV_to_J(beam->photon_energy),
		                            max_res, phases, input_items);

400
	} else {
401

Thomas White's avatar
Thomas White committed
402
		ideal_ref = new_list_intensity();
403
		input_items = read_reflections(input, ideal_ref, phases,
404
		                               NULL, esds);
405
		free(input);
Thomas White's avatar
Thomas White committed
406
407
408
409
410
		if ( check_symmetry(input_items, mero) ) {
			ERROR("The input reflection list does not appear to"
			      " have symmetry %s\n", mero);
			return 1;
		}
411

412
	}
413

414
415
	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
416

417
418
419
	if ( holo != NULL ) {
		ReflItemList *new;
		STATUS("Twinning from %s into %s\n", mero, holo);
420
421
		new = twin_reflections(ideal_ref, input_items,
		                       holo, mero, esds);
422
423
		delete_items(input_items);
		input_items = new;
Thomas White's avatar
Thomas White committed
424
425
	}

Thomas White's avatar
Thomas White committed
426
427
428
429
430
431
432
433
	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;
	}

434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
	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
453
454
455
456
	if ( template ) {
		/* Write out only reflections which are in the template
		 * (and which we have in the input) */
		ReflItemList *template_items;
457
458
		template_items = read_reflections(template,
		                                  NULL, NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
459
460
461
462
463
		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
464
		/* (quick way of copying a list) */
Thomas White's avatar
Thomas White committed
465
466
467
		union_items(write_items, input_items);
	}

468
	write_reflections(output, write_items, ideal_ref, esds, phases,
469
	                  NULL, cell);
Thomas White's avatar
Thomas White committed
470
471
472

	delete_items(input_items);
	delete_items(write_items);
473
474
475

	return 0;
}