get_hkl.c 7.59 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
29
30
31
32
33
34
35


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
"Write idealised intensity lists.\n"
"\n"
Thomas White's avatar
Thomas White committed
36
37
"  -h, --help                 Display this help message.\n"
"\n"
Thomas White's avatar
Thomas White committed
38
39
"  -t, --template=<filename>  Only include reflections mentioned in file.\n"
"      --poisson              Simulate Poisson samples.\n"
40
41
42
"  -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"
43
44
"  -o, --output=<filename>    Output filename (default: stdout).\n"
"  -i, --intensities=<file>   Read intensities from file instead of\n"
45
46
47
"                              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"
Thomas White's avatar
Thomas White committed
48
);
Thomas White's avatar
Thomas White committed
49
50
51
52
53
54
55
56
57
58
59
60
61
}


/* Apply Poisson noise to all reflections */
static void noisify_reflections(double *ref)
{
	signed int h, k, l;

	for ( h=-INDMAX; h<INDMAX; h++ ) {
	for ( k=-INDMAX; k<INDMAX; k++ ) {
	for ( l=-INDMAX; l<INDMAX; l++ ) {

		double val;
Thomas White's avatar
Thomas White committed
62
		int c;
Thomas White's avatar
Thomas White committed
63
64

		val = lookup_intensity(ref, h, k, l);
Thomas White's avatar
Thomas White committed
65
66
		c = poisson_noise(val);
		set_intensity(ref, h, k, l, c);
Thomas White's avatar
Thomas White committed
67
68
69

	}
	}
Thomas White's avatar
Thomas White committed
70
	progress_bar(h+INDMAX, 2*INDMAX, "Simulating noise");
Thomas White's avatar
Thomas White committed
71
	}
72
73
74
}


75
76
77
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
106
107
108
109
110
111
112
113
114
115
116
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
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
188
189
190
191
192
193
194
195
196
197
198
199
200
static void scold_user_about_symmetry(signed int h, signed int k, signed int l,
                                      signed int he, signed int ke,
                                      signed int le)
{
	ERROR("Merohedrally equivalent reflection (%i %i %i) found for "
	      "%i %i %i.\n", he, ke, le, h, k, l);
	ERROR("This indicates that you lied to me about the symmetry of the "
	      "input reflections.  ");
	ERROR("I won't be able to give you a meaningful result in this "
	      "situation, so I'm going to give up right now.  ");
	ERROR("Please reconsider your previous processing of the data, and "
	      "perhaps try again with a lower symmetry for the '-y' option.\n");
}


static int find_unique_equiv(ReflItemList *items, signed int h, signed int k,
                             signed int l, const char *mero, signed int *hu,
                             signed int *ku, signed int *lu)
{
	int i;

	for ( i=0; i<num_equivs(h, k, l, mero); i++ ) {

		signed int he, ke, le;
		get_equiv(h, k, l, &he, &ke, &le, mero, i);
		if ( find_item(items, he, ke, le) ) {
			*hu = he;  *ku = ke;  *lu = le;
			return 1;
		}

	}

	return 0;
}


static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
                                      const char *holo, const char *mero)
{
	int i;
	ReflItemList *new;

	new = new_items();

	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);
		h = it->h;  k = it->k;  l = it->l;

		/* None of the equivalent reflections should exist in the
		 * input dataset.  That would indicate that the user lied about
		 * the input symmetry.
		 *
		 * Start from j=1 to ignore the reflection itself.
		 */
		for ( j=1; j<num_equivs(h, k, l, mero); j++ ) {

			signed int he, ke, le;
			get_equiv(h, k, l, &he, &ke, &le, mero, j);
			if ( !find_item(items, he, ke, le) ) continue;

			scold_user_about_symmetry(h, k, l, he, ke, le);
			abort();

		}
		/* It doesn't matter if the reflection wasn't actually the one
		 * we define as being in the asymmetric unit cell, as long as
		 * things aren't confused by there being more than one of it.
		 */

		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.
			 *
			 * We checked earlier that there's only one of these
			 * for each reflection.
			 */
			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;
}


201
202
203
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
204
	double *ideal_ref;
205
	double *phases;
206
	struct molecule *mol;
Thomas White's avatar
Thomas White committed
207
208
	char *template = NULL;
	int config_noisify = 0;
209
210
	char *holo = NULL;
	char *mero = NULL;
Thomas White's avatar
Thomas White committed
211
	char *output = NULL;
212
	char *input = NULL;
213
	char *filename = NULL;
Thomas White's avatar
Thomas White committed
214
215
	ReflItemList *input_items;
	ReflItemList *write_items;
216
217
218
219

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
220
221
		{"template",           1, NULL,               't'},
		{"poisson",            0, &config_noisify,     1},
Thomas White's avatar
Thomas White committed
222
		{"output",             1, NULL,               'o'},
223
224
		{"twin",               1, NULL,               'w'},
		{"symmetry",           1, NULL,               'y'},
225
		{"intensities",        1, NULL,               'i'},
226
		{"pdb",                1, NULL,               'p'},
227
228
229
230
		{0, 0, NULL, 0}
	};

	/* Short options */
231
	while ((c = getopt_long(argc, argv, "ht:o:i:p:w:y:", longopts, NULL)) != -1) {
232
233

		switch (c) {
Thomas White's avatar
Thomas White committed
234
		case 'h' :
235
236
237
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
238
		case 't' :
Thomas White's avatar
Thomas White committed
239
240
241
			template = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
242
		case 'o' :
Thomas White's avatar
Thomas White committed
243
244
245
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
246
		case 'i' :
247
248
249
			input = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
250
		case 'p' :
251
252
253
			filename = strdup(optarg);
			break;

254
255
256
257
258
259
260
261
		case 'y' :
			mero = strdup(optarg);
			break;

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

Thomas White's avatar
Thomas White committed
262
		case 0 :
263
264
			break;

Thomas White's avatar
Thomas White committed
265
		default :
266
267
268
269
270
			return 1;
		}

	}

271
272
273
274
275
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

	mol = load_molecule(filename);
Thomas White's avatar
Thomas White committed
276
	phases = new_list_phase();
277
	if ( input == NULL ) {
Thomas White's avatar
Thomas White committed
278
		input_items = new_items();
279
		ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9),
Thomas White's avatar
Thomas White committed
280
		                            phases, input_items);
281
	} else {
Thomas White's avatar
Thomas White committed
282
283
284
		ideal_ref = new_list_intensity();
		phases = new_list_phase();
		input_items = read_reflections(input, ideal_ref, phases, NULL);
285
286
		free(input);
	}
287

288
	if ( config_noisify ) noisify_reflections(ideal_ref);
Thomas White's avatar
Thomas White committed
289

290
291
292
293
294
295
	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
296
297
	}

Thomas White's avatar
Thomas White committed
298
299
300
301
302
303
304
305
306
307
	if ( template ) {
		/* Write out only reflections which are in the template
		 * (and which we have in the input) */
		ReflItemList *template_items;
		template_items = read_reflections(template, NULL, NULL, NULL);
		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
308
		/* (quick way of copying a list) */
Thomas White's avatar
Thomas White committed
309
310
311
312
313
314
315
316
		union_items(write_items, input_items);
	}

	write_reflections(output, write_items, ideal_ref, phases, NULL,
	                  mol->cell);

	delete_items(input_items);
	delete_items(write_items);
317
318
319

	return 0;
}