get_hkl.c 6.03 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
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);

93
94
95
96
		/* 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.
97
		 */
98
99
		get_asymm(it->h, it->k, it->l, &h, &k, &l, holo);
		if ( find_item(new, h, k, l) ) continue;
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

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


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

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
167
168
		{"template",           1, NULL,               't'},
		{"poisson",            0, &config_noisify,     1},
Thomas White's avatar
Thomas White committed
169
		{"output",             1, NULL,               'o'},
170
171
		{"twin",               1, NULL,               'w'},
		{"symmetry",           1, NULL,               'y'},
172
		{"intensities",        1, NULL,               'i'},
173
		{"pdb",                1, NULL,               'p'},
174
175
176
177
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
181
		case 'h' :
182
183
184
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
185
		case 't' :
Thomas White's avatar
Thomas White committed
186
187
188
			template = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
189
		case 'o' :
Thomas White's avatar
Thomas White committed
190
191
192
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
193
		case 'i' :
194
195
196
			input = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
197
		case 'p' :
198
199
200
			filename = strdup(optarg);
			break;

201
202
203
204
205
206
207
208
		case 'y' :
			mero = strdup(optarg);
			break;

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

Thomas White's avatar
Thomas White committed
209
		case 0 :
210
211
			break;

Thomas White's avatar
Thomas White committed
212
		default :
213
214
215
216
217
			return 1;
		}

	}

218
219
220
221
222
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

	mol = load_molecule(filename);
Thomas White's avatar
Thomas White committed
223
	phases = new_list_phase();
224
	if ( input == NULL ) {
Thomas White's avatar
Thomas White committed
225
		input_items = new_items();
226
		ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9),
Thomas White's avatar
Thomas White committed
227
		                            phases, input_items);
228
	} else {
Thomas White's avatar
Thomas White committed
229
230
231
		ideal_ref = new_list_intensity();
		phases = new_list_phase();
		input_items = read_reflections(input, ideal_ref, phases, NULL);
232
233
		free(input);
	}
234

235
	if ( config_noisify ) noisify_reflections(ideal_ref);
Thomas White's avatar
Thomas White committed
236

237
238
239
240
241
242
	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
243
244
	}

Thomas White's avatar
Thomas White committed
245
246
247
248
249
250
251
252
253
254
	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
255
		/* (quick way of copying a list) */
Thomas White's avatar
Thomas White committed
256
257
258
259
260
261
262
263
		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);
264
265
266

	return 0;
}