get_hkl.c 6.9 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"
48
49
50
"      --no-phases            Do not try to use phases in the input file.\n"
"      --multiplicity         Multiply intensities by the number of\n"
"                              equivalent reflections.\n"
Thomas White's avatar
Thomas White committed
51
);
Thomas White's avatar
Thomas White committed
52
53
54
55
56
57
58
59
60
61
62
63
64
}


/* 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
65
		int c;
Thomas White's avatar
Thomas White committed
66
67

		val = lookup_intensity(ref, h, k, l);
Thomas White's avatar
Thomas White committed
68
69
		c = poisson_noise(val);
		set_intensity(ref, h, k, l, c);
Thomas White's avatar
Thomas White committed
70
71
72

	}
	}
Thomas White's avatar
Thomas White committed
73
	progress_bar(h+INDMAX, 2*INDMAX, "Simulating noise");
Thomas White's avatar
Thomas White committed
74
	}
75
76
77
}


78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
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);

96
97
98
99
		/* There is a many-to-one correspondence between reflections
		 * in the merohedral and holohedral unit cells.  Do the
		 * calculation only once for each reflection in the holohedral
		 * cell, which contains fewer reflections.
100
		 */
101
102
		get_asymm(it->h, it->k, it->l, &h, &k, &l, holo);
		if ( find_item(new, h, k, l) ) continue;
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

		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;
}


151
152
153
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
154
	double *ideal_ref;
155
	double *phases;
156
	struct molecule *mol;
Thomas White's avatar
Thomas White committed
157
158
	char *template = NULL;
	int config_noisify = 0;
159
160
	int config_nophase = 0;
	int config_multi = 0;
161
162
	char *holo = NULL;
	char *mero = NULL;
Thomas White's avatar
Thomas White committed
163
	char *output = NULL;
164
	char *input = NULL;
165
	char *filename = NULL;
Thomas White's avatar
Thomas White committed
166
167
	ReflItemList *input_items;
	ReflItemList *write_items;
168
169
170
171

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
172
173
		{"template",           1, NULL,               't'},
		{"poisson",            0, &config_noisify,     1},
Thomas White's avatar
Thomas White committed
174
		{"output",             1, NULL,               'o'},
175
176
		{"twin",               1, NULL,               'w'},
		{"symmetry",           1, NULL,               'y'},
177
		{"intensities",        1, NULL,               'i'},
178
		{"pdb",                1, NULL,               'p'},
179
180
		{"no-phases",          0, &config_nophase,     1},
		{"multiplicity",       0, &config_multi,       1},
181
182
183
184
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
188
		case 'h' :
189
190
191
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
192
		case 't' :
Thomas White's avatar
Thomas White committed
193
194
195
			template = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
196
		case 'o' :
Thomas White's avatar
Thomas White committed
197
198
199
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
200
		case 'i' :
201
202
203
			input = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
204
		case 'p' :
205
206
207
			filename = strdup(optarg);
			break;

208
209
210
211
212
213
214
215
		case 'y' :
			mero = strdup(optarg);
			break;

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

Thomas White's avatar
Thomas White committed
216
		case 0 :
217
218
			break;

Thomas White's avatar
Thomas White committed
219
		default :
220
221
222
223
224
			return 1;
		}

	}

225
226
227
228
229
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

	mol = load_molecule(filename);
230
231
232
233
234
	if ( !config_nophase ) {
		phases = new_list_phase();
	} else {
		phases = NULL;
	}
235
	if ( input == NULL ) {
Thomas White's avatar
Thomas White committed
236
		input_items = new_items();
237
		ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9),
Thomas White's avatar
Thomas White committed
238
		                            phases, input_items);
239
	} else {
Thomas White's avatar
Thomas White committed
240
		ideal_ref = new_list_intensity();
241
242
		input_items = read_reflections(input, ideal_ref, phases,
		                               NULL, NULL);
243
244
		free(input);
	}
245

246
	if ( config_noisify ) noisify_reflections(ideal_ref);
Thomas White's avatar
Thomas White committed
247

248
249
250
251
252
253
	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
254
255
	}

256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
	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
275
276
277
278
	if ( template ) {
		/* Write out only reflections which are in the template
		 * (and which we have in the input) */
		ReflItemList *template_items;
279
280
		template_items = read_reflections(template,
		                                  NULL, NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
281
282
283
284
285
		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
286
		/* (quick way of copying a list) */
Thomas White's avatar
Thomas White committed
287
288
289
290
291
292
293
294
		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);
295
296
297

	return 0;
}