1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.crystal;
39
40 import org.apache.commons.configuration2.CompositeConfiguration;
41
42 import javax.annotation.Nullable;
43 import java.util.ArrayList;
44 import java.util.HashMap;
45 import java.util.Map.Entry;
46 import java.util.logging.Logger;
47
48 import static ffx.crystal.SymOp.applySymRot;
49 import static ffx.crystal.SymOp.applyTransSymRot;
50 import static ffx.numerics.math.ScalarMath.mod;
51 import static org.apache.commons.math3.util.FastMath.PI;
52 import static org.apache.commons.math3.util.FastMath.cos;
53 import static org.apache.commons.math3.util.FastMath.floor;
54 import static org.apache.commons.math3.util.FastMath.max;
55 import static org.apache.commons.math3.util.FastMath.min;
56 import static org.apache.commons.math3.util.FastMath.rint;
57
58
59
60
61
62
63
64
65
66
67 public class ReflectionList {
68
69 private static final Logger logger = Logger.getLogger(ReflectionList.class.getName());
70
71
72
73
74 public final ArrayList<HKL> hklList = new ArrayList<>();
75
76
77
78 public final Crystal crystal;
79
80
81
82 public final SpaceGroup spaceGroup;
83
84
85
86 public final Resolution resolution;
87
88
89
90 final HashMap<String, HKL> hklMap = new HashMap<>();
91
92
93
94 public int nBins = 10;
95
96
97
98 private final double[] histogram = new double[1001];
99
100
101
102 private double minResolution;
103
104
105
106 private double maxResolution;
107
108
109
110
111
112
113
114 public ReflectionList(Crystal crystal, Resolution resolution) {
115 this(crystal, resolution, null);
116 }
117
118
119
120
121
122
123
124
125 public ReflectionList(Crystal crystal, Resolution resolution, CompositeConfiguration properties) {
126 this.crystal = crystal;
127 this.spaceGroup = crystal.spaceGroup;
128 this.resolution = resolution;
129
130 int hMax = (int) (this.crystal.a / this.resolution.resolutionLimit());
131 int kMax = (int) (this.crystal.b / this.resolution.resolutionLimit());
132 int lMax = (int) (this.crystal.c / this.resolution.resolutionLimit());
133
134 minResolution = Double.POSITIVE_INFINITY;
135 maxResolution = Double.NEGATIVE_INFINITY;
136 int n = 0;
137
138 HKL hkl = new HKL();
139 for (int h = -hMax; h <= hMax; h++) {
140 hkl.setH(h);
141 for (int k = -kMax; k <= kMax; k++) {
142 hkl.setK(k);
143 for (int l = -lMax; l <= lMax; l++) {
144 hkl.setL(l);
145
146 double res = this.crystal.invressq(hkl);
147 getEpsilon(hkl);
148 LaueSystem laueSystem = spaceGroup.laueSystem;
149 if (laueSystem.checkRestrictions(h, k, l)
150 && resolution.inInverseResSqRange(res)
151 && !hkl.sysAbs()) {
152 minResolution = min(res, minResolution);
153 maxResolution = max(res, maxResolution);
154 String s = (h + "_" + k + "_" + l).intern();
155 hklMap.put(s, new HKL(hkl.getH(), hkl.getK(), hkl.getL(), hkl.getEpsilon(), hkl.allowed));
156 n++;
157 }
158 }
159 }
160 }
161
162 n = 0;
163 for (Entry<String, HKL> entry : hklMap.entrySet()) {
164 HKL ih = entry.getValue();
165 ih.setIndex(n);
166 hklList.add(ih);
167 n++;
168 }
169
170
171 for (HKL ih : hklList) {
172 double r = (this.crystal.invressq(ih) - minResolution) / (maxResolution - minResolution);
173 int i = (int) (min(r, 0.999) * 1000.0);
174 histogram[i + 1] += 1.0;
175 }
176
177
178 for (int i = 1; i < histogram.length; i++) {
179 histogram[i] += histogram[i - 1];
180 }
181 for (int i = 0; i < histogram.length; i++) {
182 histogram[i] /= histogram[histogram.length - 1];
183 }
184
185
186 setResolutionBins(properties);
187 }
188
189
190
191
192
193
194
195
196
197
198
199
200
201 public ReflectionList(double a, double b, double c, double alpha, double beta, double gamma,
202 String sg, double resolution) {
203 this(new Crystal(a, b, c, alpha, beta, gamma, sg), new Resolution(resolution));
204 }
205
206
207
208
209
210
211
212
213 public boolean findSymHKL(HKL hkl, HKL mate) {
214 return findSymHKL(hkl, mate, false);
215 }
216
217
218
219
220
221
222
223
224
225
226 public boolean findSymHKL(int h, int k, int l, HKL mate) {
227 return findSymHKL(new HKL(h, k, l), mate, false);
228 }
229
230
231
232
233
234
235
236
237
238
239
240 public boolean findSymHKL(int h, int k, int l, HKL mate, boolean transpose) {
241 return findSymHKL(new HKL(h, k, l), mate, transpose);
242 }
243
244
245
246
247
248
249
250
251
252 public HKL getHKL(int h, int k, int l) {
253 String s = (h + "_" + k + "_" + l);
254 return hklMap.get(s);
255 }
256
257
258
259
260
261
262
263 public HKL getHKL(HKL hkl) {
264 return getHKL(hkl.getH(), hkl.getK(), hkl.getL());
265 }
266
267
268
269
270
271
272 public double getMaxResolution() {
273 return maxResolution;
274 }
275
276
277
278
279
280
281 public double getMinResolution() {
282 return minResolution;
283 }
284
285
286
287
288
289
290
291 public boolean hasHKL(HKL hkl) {
292 return hasHKL(hkl.getH(), hkl.getK(), hkl.getL());
293 }
294
295
296
297
298
299
300
301 public final double ordinal(double s) {
302 double r = (s - minResolution) / (maxResolution - minResolution);
303 r = min(r, 0.999) * 1000.0;
304 int i = (int) r;
305 r -= floor(r);
306 return ((1.0 - r) * histogram[i] + r * histogram[i + 1]);
307 }
308
309
310
311
312 @Override
313 public String toString() {
314 return " Reflection list with "
315 + this.hklList.size()
316 + " reflections, spacegroup "
317 + this.spaceGroup.shortName
318 + " resolution limit: "
319 + resolution.resolutionLimit();
320 }
321
322
323
324
325
326
327
328
329
330 private boolean findSymHKL(HKL hkl, HKL mate, boolean transpose) {
331 int nsym = spaceGroup.numPrimitiveSymEquiv;
332
333 for (int i = 0; i < nsym; i++) {
334 if (transpose) {
335 applyTransSymRot(hkl, mate, spaceGroup.symOps.get(i));
336 } else {
337 applySymRot(hkl, mate, spaceGroup.symOps.get(i));
338 }
339
340 LaueSystem laueSystem = spaceGroup.laueSystem;
341 if (laueSystem.checkRestrictions(mate.getH(), mate.getK(), mate.getL())) {
342 return false;
343 }
344 if (laueSystem.checkRestrictions(-mate.getH(), -mate.getK(), -mate.getL())) {
345 mate.setH(-mate.getH());
346 mate.setK(-mate.getK());
347 mate.setL(-mate.getL());
348 return true;
349 }
350 }
351
352 mate.setH(hkl.getH());
353 mate.setK(hkl.getK());
354 mate.setL(hkl.getL());
355 return false;
356 }
357
358
359
360
361
362
363
364
365
366 private boolean hasHKL(int h, int k, int l) {
367 String s = (h + "_" + k + "_" + l);
368 return hklMap.containsKey(s);
369 }
370
371 private void getEpsilon(HKL hkl) {
372 int epsilon = 1;
373 int allowed = 255;
374
375 int nSym = spaceGroup.symOps.size();
376 for (int i = 1; i < nSym; i++) {
377 HKL mate = new HKL();
378 SymOp symOp = spaceGroup.symOps.get(i);
379 applySymRot(hkl, mate, symOp);
380 double shift = symOp.symPhaseShift(hkl);
381 if (mate.equals(hkl)) {
382 if (cos(shift) > 0.999) {
383 epsilon++;
384 } else {
385 allowed = 0;
386 epsilon = 0;
387 break;
388 }
389 } else if (mate.equals(hkl.neg())) {
390
391 allowed = (int) rint(mod(-0.5 * shift, PI) / (PI / HKL.ndiv));
392 }
393 }
394 if (hkl.getH() == 0 && hkl.getK() == 0 && hkl.getL() == 0) {
395 allowed = 0;
396 }
397
398 hkl.setEpsilon(epsilon);
399 hkl.setAllowed(allowed);
400 }
401
402 private void setResolutionBins(@Nullable CompositeConfiguration properties) {
403 if (properties != null) {
404 nBins = properties.getInt("reflection-bins", 10);
405 }
406 double nBinsDouble = nBins;
407 for (HKL ih : hklList) {
408 int bin = (int) (nBinsDouble * ordinal(crystal.invressq(ih)));
409 ih.bin = min(bin, nBins - 1);
410 }
411 }
412 }