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.xray.parsers;
39
40 import ffx.crystal.Crystal;
41 import ffx.crystal.HKL;
42 import ffx.crystal.ReflectionList;
43 import ffx.crystal.Resolution;
44 import ffx.crystal.SpaceGroupInfo;
45 import ffx.numerics.math.ComplexNumber;
46 import ffx.xray.DiffractionRefinementData;
47 import ffx.xray.parsers.MTZWriter.MTZType;
48 import org.apache.commons.configuration2.CompositeConfiguration;
49 import org.apache.commons.lang3.Strings;
50
51 import java.io.DataInputStream;
52 import java.io.EOFException;
53 import java.io.File;
54 import java.io.FileInputStream;
55 import java.io.IOException;
56 import java.nio.ByteBuffer;
57 import java.nio.ByteOrder;
58 import java.util.ArrayList;
59 import java.util.Iterator;
60 import java.util.logging.Level;
61 import java.util.logging.Logger;
62
63 import static java.lang.Double.parseDouble;
64 import static java.lang.Float.parseFloat;
65 import static java.lang.Integer.parseInt;
66 import static java.lang.String.format;
67 import static org.apache.commons.math3.util.FastMath.cos;
68 import static org.apache.commons.math3.util.FastMath.max;
69 import static org.apache.commons.math3.util.FastMath.min;
70 import static org.apache.commons.math3.util.FastMath.sin;
71 import static org.apache.commons.math3.util.FastMath.sqrt;
72 import static org.apache.commons.math3.util.FastMath.toRadians;
73
74
75
76
77
78
79
80
81
82
83 public class MTZFilter implements DiffractionFileFilter {
84
85 private static final Logger logger = Logger.getLogger(MTZFilter.class.getName());
86 private final ArrayList<Column> columns = new ArrayList<>();
87 private final ArrayList<Dataset> dataSets = new ArrayList<>();
88 private boolean headerParsed = false;
89 private String title;
90 private String foString, sigFoString, rFreeString;
91 private int h, k, l, fo, sigFo, rFree;
92 private int fPlus, sigFPlus, fMinus, sigFMinus, rFreePlus, rFreeMinus;
93 private int fc, phiC, fs, phiS;
94 private int dsetOffset = 1;
95 private int nColumns;
96 private int nReflections;
97 private int spaceGroupNum;
98 private String spaceGroupName;
99 private double resLow;
100 private double resHigh;
101
102
103 public MTZFilter() {
104 }
105
106
107
108
109
110
111
112
113
114
115 public void averageFcs(File mtzFile1, File mtzFile2, ReflectionList reflectionlist,
116 int iter, CompositeConfiguration properties) {
117
118 DiffractionRefinementData fcdata1 = new DiffractionRefinementData(properties, reflectionlist);
119 DiffractionRefinementData fcdata2 = new DiffractionRefinementData(properties, reflectionlist);
120
121 readFcs(mtzFile1, reflectionlist, fcdata1, properties);
122 readFcs(mtzFile2, reflectionlist, fcdata2, properties);
123
124
125 logger.info(format(" Iteration for averaging: %d.", iter));
126 for (int i = 0; i < reflectionlist.hklList.size(); i++) {
127 fcdata1.fc[i][0] += (fcdata2.fc[i][0] - fcdata1.fc[i][0]) / iter;
128 fcdata1.fc[i][1] += (fcdata2.fc[i][1] - fcdata1.fc[i][1]) / iter;
129
130 fcdata1.fs[i][0] += (fcdata2.fs[i][0] - fcdata1.fs[i][0]) / iter;
131 fcdata1.fs[i][1] += (fcdata2.fs[i][1] - fcdata1.fs[i][1]) / iter;
132 }
133
134
135 MTZWriter mtzOut = new MTZWriter(reflectionlist, fcdata1, mtzFile1.getName(), MTZType.FCONLY);
136 mtzOut.write();
137 }
138
139
140 @Override
141 public ReflectionList getReflectionList(File mtzFile, CompositeConfiguration properties) {
142 ByteOrder byteOrder = ByteOrder.nativeOrder();
143 FileInputStream fileInputStream;
144 DataInputStream dataInputStream;
145 try {
146 fileInputStream = new FileInputStream(mtzFile);
147 dataInputStream = new DataInputStream(fileInputStream);
148
149 byte[] headerOffset = new byte[4];
150 byte[] bytes = new byte[80];
151 int offset = 0;
152
153
154 dataInputStream.read(bytes, offset, 4);
155
156
157 dataInputStream.read(headerOffset, offset, 4);
158
159
160 dataInputStream.read(bytes, offset, 4);
161 ByteBuffer byteBuffer = ByteBuffer.wrap(bytes);
162 int stamp = byteBuffer.order(ByteOrder.BIG_ENDIAN).getInt();
163 String stampstr = Integer.toHexString(stamp);
164 switch (stampstr.charAt(0)) {
165 case '1':
166 case '3':
167 if (byteOrder.equals(ByteOrder.LITTLE_ENDIAN)) {
168 byteOrder = ByteOrder.BIG_ENDIAN;
169 }
170 break;
171 case '4':
172 if (byteOrder.equals(ByteOrder.BIG_ENDIAN)) {
173 byteOrder = ByteOrder.LITTLE_ENDIAN;
174 }
175 break;
176 }
177
178 byteBuffer = ByteBuffer.wrap(headerOffset);
179 int headerOffsetI = byteBuffer.order(byteOrder).getInt();
180
181
182 dataInputStream.skipBytes((headerOffsetI - 4) * 4);
183
184 for (Boolean parsing = true; parsing; dataInputStream.read(bytes, offset, 80)) {
185 String mtzstr = new String(bytes);
186 parsing = parseHeader(mtzstr);
187 }
188 } catch (EOFException e) {
189 String message = " MTZ end of file reached.";
190 logger.log(Level.WARNING, message, e);
191 return null;
192 } catch (IOException e) {
193 String message = " MTZ IO exception.";
194 logger.log(Level.WARNING, message, e);
195 return null;
196 }
197
198
199 foString = sigFoString = rFreeString = null;
200 if (properties != null) {
201 foString = properties.getString("fostring", null);
202 sigFoString = properties.getString("sigfostring", null);
203 rFreeString = properties.getString("rfreestring", null);
204 }
205 h = k = l = fo = sigFo = rFree = -1;
206 fPlus = sigFPlus = fMinus = sigFMinus = rFreePlus = rFreeMinus = -1;
207 fc = phiC = -1;
208 boolean print = false;
209 parseColumns(print);
210 parseFcColumns(print);
211
212 if (fo < 0 && fPlus < 0 && sigFo < 0 && sigFPlus < 0 && fc < 0 && phiC < 0) {
213 logger.info(
214 " The MTZ header contains insufficient information to generate the reflection list.");
215 logger.info(
216 " For non-default column labels set fostring/sigfostring in the properties file.");
217 return null;
218 }
219
220 Column column;
221 if (fo > 0) {
222 column = columns.get(fo);
223 } else if (fPlus > 0) {
224 column = columns.get(fPlus);
225 } else {
226 column = columns.get(fc);
227 }
228 Dataset dataSet = dataSets.get(column.id - dsetOffset);
229
230 if (logger.isLoggable(Level.INFO)) {
231 StringBuilder sb = new StringBuilder();
232 sb.append(format(" Reading %s\n", mtzFile.getName()));
233 sb.append(" Setting up reflection list based on MTZ file.\n");
234 sb.append(format(" Space group number: %d (name: %s)\n",
235 spaceGroupNum, SpaceGroupInfo.spaceGroupNames[spaceGroupNum - 1]));
236 sb.append(format(" Resolution: %8.3f\n", 0.999999 * resHigh));
237 sb.append(format(" Cell: %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f",
238 dataSet.cell[0], dataSet.cell[1], dataSet.cell[2],
239 dataSet.cell[3], dataSet.cell[4], dataSet.cell[5]));
240 logger.info(sb.toString());
241 }
242
243 Crystal crystal = new Crystal(dataSet.cell[0], dataSet.cell[1], dataSet.cell[2],
244 dataSet.cell[3], dataSet.cell[4], dataSet.cell[5],
245 SpaceGroupInfo.spaceGroupNames[spaceGroupNum - 1]);
246
247 double sampling = 0.6;
248 if (properties != null) {
249 sampling = properties.getDouble("sampling", 0.6);
250 }
251 Resolution resolution = new Resolution(0.999999 * resHigh, sampling);
252 return new ReflectionList(crystal, resolution, properties);
253 }
254
255
256 @Override
257 public ReflectionList getReflectionList(File mtzFile) {
258 return getReflectionList(mtzFile, null);
259 }
260
261
262 @Override
263 public double getResolution(File mtzFile, Crystal crystal) {
264 getReflectionList(mtzFile, null);
265 return resHigh;
266 }
267
268
269 public void printHeader() {
270 if (logger.isLoggable(Level.INFO)) {
271 StringBuilder sb = new StringBuilder();
272 sb.append(" MTZ title: ").append(title).append("\n");
273 sb.append(" MTZ space group: ")
274 .append(spaceGroupName)
275 .append(" space group number: ")
276 .append(spaceGroupNum)
277 .append(" (")
278 .append(SpaceGroupInfo.spaceGroupNames[spaceGroupNum - 1])
279 .append(")\n");
280 sb.append(" MTZ resolution: ").append(resLow).append(" - ").append(resHigh).append("\n");
281 sb.append(" Number of reflections: ").append(nReflections).append("\n");
282
283 int ndset = 1;
284 for (Iterator<Dataset> i = dataSets.iterator(); i.hasNext(); ndset++) {
285 Dataset d = i.next();
286 sb.append(" dataset ").append(ndset).append(": ").append(d.dataset).append("\n");
287 sb.append(" project ").append(ndset).append(": ").append(d.project).append("\n");
288 sb.append(" wavelength ").append(ndset).append(": ").append(d.lambda).append("\n");
289 sb.append(" cell ").append(ndset).append(": ")
290 .append(d.cell[0]).append(" ").append(d.cell[1]).append(" ")
291 .append(d.cell[2]).append(" ").append(d.cell[3]).append(" ")
292 .append(d.cell[4]).append(" ").append(d.cell[5]).append("\n");
293 sb.append("\n");
294 }
295
296 sb.append(" Number of columns: ").append(nColumns).append("\n");
297 int nc = 0;
298 for (Iterator<Column> i = columns.iterator(); i.hasNext(); nc++) {
299 Column c = i.next();
300 sb.append(format(" column %d: dataset id: %d min: %9.2f max: %9.2f label: %s type: %c\n",
301 nc, c.id, c.min, c.max, c.label, c.type));
302 }
303 logger.info(sb.toString());
304 }
305 }
306
307
308 @Override
309 public boolean readFile(File mtzFile, ReflectionList reflectionList,
310 DiffractionRefinementData refinementData, CompositeConfiguration properties) {
311 int nRead, nIgnore, nRes, nFriedel, nCut;
312 ByteOrder byteOrder = ByteOrder.nativeOrder();
313 FileInputStream fileInputStream;
314 DataInputStream dataInputStream;
315 boolean transpose = false;
316
317 StringBuilder sb = new StringBuilder();
318 try {
319 fileInputStream = new FileInputStream(mtzFile);
320 dataInputStream = new DataInputStream(fileInputStream);
321
322 byte[] headerOffset = new byte[4];
323 byte[] bytes = new byte[80];
324 int offset = 0;
325
326
327 dataInputStream.read(bytes, offset, 4);
328
329
330 dataInputStream.read(headerOffset, offset, 4);
331
332
333 dataInputStream.read(bytes, offset, 4);
334 ByteBuffer byteBuffer = ByteBuffer.wrap(bytes);
335 int stamp = byteBuffer.order(ByteOrder.BIG_ENDIAN).getInt();
336 String stampString = Integer.toHexString(stamp);
337 switch (stampString.charAt(0)) {
338 case '1':
339 case '3':
340 if (byteOrder.equals(ByteOrder.LITTLE_ENDIAN)) {
341 byteOrder = ByteOrder.BIG_ENDIAN;
342 }
343 break;
344 case '4':
345 if (byteOrder.equals(ByteOrder.BIG_ENDIAN)) {
346 byteOrder = ByteOrder.LITTLE_ENDIAN;
347 }
348 break;
349 }
350
351 byteBuffer = ByteBuffer.wrap(headerOffset);
352 int headerOffsetI = byteBuffer.order(byteOrder).getInt();
353
354
355 dataInputStream.skipBytes((headerOffsetI - 4) * 4);
356
357 for (Boolean parsing = true; parsing; dataInputStream.read(bytes, offset, 80)) {
358 String mtzstr = new String(bytes);
359 parsing = parseHeader(mtzstr);
360 }
361
362
363 foString = sigFoString = rFreeString = null;
364 if (properties != null) {
365 foString = properties.getString("fostring", null);
366 sigFoString = properties.getString("sigfostring", null);
367 rFreeString = properties.getString("rfreestring", null);
368 }
369 h = k = l = fo = sigFo = rFree = -1;
370 fPlus = sigFPlus = fMinus = sigFMinus = rFreePlus = rFreeMinus = -1;
371 boolean print = true;
372 parseColumns(print);
373
374 if (h < 0 || k < 0 || l < 0) {
375 String message = "Fatal error in MTZ file - no H K L indexes?\n";
376 logger.log(Level.SEVERE, message);
377 return false;
378 }
379
380
381 fileInputStream = new FileInputStream(mtzFile);
382 dataInputStream = new DataInputStream(fileInputStream);
383
384
385 dataInputStream.skipBytes(80);
386
387
388 float[] data = new float[nColumns];
389 HKL mate = new HKL();
390 int nPosIgnore = 0;
391 int nTransIgnore = 0;
392 int nZero = 0;
393 int none = 0;
394 for (int i = 0; i < nReflections; i++) {
395 for (int j = 0; j < nColumns; j++) {
396 dataInputStream.read(bytes, offset, 4);
397 byteBuffer = ByteBuffer.wrap(bytes);
398 data[j] = byteBuffer.order(byteOrder).getFloat();
399 }
400 int ih = (int) data[h];
401 int ik = (int) data[k];
402 int il = (int) data[l];
403 reflectionList.findSymHKL(ih, ik, il, mate, false);
404 HKL hklpos = reflectionList.getHKL(mate);
405 if (hklpos == null) {
406 nPosIgnore++;
407 }
408
409 reflectionList.findSymHKL(ih, ik, il, mate, true);
410 HKL hkltrans = reflectionList.getHKL(mate);
411 if (hkltrans == null) {
412 nTransIgnore++;
413 }
414 if (rFree > 0) {
415 if (((int) data[rFree]) == 0) {
416 nZero++;
417 } else if (((int) data[rFree]) == 1) {
418 none++;
419 }
420 }
421 if (rFreePlus > 0) {
422 if (((int) data[rFreePlus]) == 0) {
423 nZero++;
424 } else if (((int) data[rFreePlus]) == 1) {
425 none++;
426 }
427 }
428 if (rFreeMinus > 0) {
429 if (((int) data[rFreeMinus]) == 0) {
430 nZero++;
431 } else if (((int) data[rFreeMinus]) == 1) {
432 none++;
433 }
434 }
435 }
436 if (nPosIgnore > nTransIgnore) {
437 transpose = true;
438 }
439
440 if (none > (nZero * 2) && refinementData.rFreeFlag < 0) {
441 refinementData.setFreeRFlag(0);
442 sb.append(format( " Setting R free flag to %d based on MTZ file data.\n", refinementData.rFreeFlag));
443 } else if (nZero > (none * 2) && refinementData.rFreeFlag < 0) {
444 refinementData.setFreeRFlag(1);
445 sb.append(format(" Setting R free flag to %d based on MTZ file data.\n", refinementData.rFreeFlag));
446 } else if (refinementData.rFreeFlag < 0) {
447 refinementData.setFreeRFlag(0);
448 sb.append(format(" Setting R free flag to MTZ default: %d\n", refinementData.rFreeFlag));
449 }
450
451
452 fileInputStream = new FileInputStream(mtzFile);
453 dataInputStream = new DataInputStream(fileInputStream);
454
455
456 dataInputStream.skipBytes(80);
457
458
459 double[][] anofSigF = new double[refinementData.n][4];
460 for (int i = 0; i < refinementData.n; i++) {
461 anofSigF[i][0] = anofSigF[i][1] = anofSigF[i][2] = anofSigF[i][3] = Double.NaN;
462 }
463 nRead = nIgnore = nRes = nFriedel = nCut = 0;
464 for (int i = 0; i < nReflections; i++) {
465 for (int j = 0; j < nColumns; j++) {
466 dataInputStream.read(bytes, offset, 4);
467 byteBuffer = ByteBuffer.wrap(bytes);
468 data[j] = byteBuffer.order(byteOrder).getFloat();
469 }
470 int ih = (int) data[h];
471 int ik = (int) data[k];
472 int il = (int) data[l];
473 boolean friedel = reflectionList.findSymHKL(ih, ik, il, mate, transpose);
474 HKL hkl = reflectionList.getHKL(mate);
475 if (hkl != null) {
476 if (fo > 0 && sigFo > 0) {
477 if (refinementData.fSigFCutoff > 0.0) {
478 if ((data[fo] / data[sigFo]) < refinementData.fSigFCutoff) {
479 nCut++;
480 continue;
481 }
482 }
483 if (friedel) {
484 anofSigF[hkl.getIndex()][2] = data[fo];
485 anofSigF[hkl.getIndex()][3] = data[sigFo];
486 nFriedel++;
487 } else {
488 anofSigF[hkl.getIndex()][0] = data[fo];
489 anofSigF[hkl.getIndex()][1] = data[sigFo];
490 }
491 } else {
492 if (fPlus > 0 && sigFPlus > 0) {
493 if (refinementData.fSigFCutoff > 0.0) {
494 if ((data[fPlus] / data[sigFPlus]) < refinementData.fSigFCutoff) {
495 nCut++;
496 continue;
497 }
498 }
499 anofSigF[hkl.getIndex()][0] = data[fPlus];
500 anofSigF[hkl.getIndex()][1] = data[sigFPlus];
501 }
502 if (fMinus > 0 && sigFMinus > 0) {
503 if (refinementData.fSigFCutoff > 0.0) {
504 if ((data[fMinus] / data[sigFMinus]) < refinementData.fSigFCutoff) {
505 nCut++;
506 continue;
507 }
508 }
509 anofSigF[hkl.getIndex()][2] = data[fMinus];
510 anofSigF[hkl.getIndex()][3] = data[sigFMinus];
511 }
512 }
513 if (rFree > 0) {
514 refinementData.setFreeR(hkl.getIndex(), (int) data[rFree]);
515 } else {
516 if (rFreePlus > 0 && rFreeMinus > 0) {
517
518 refinementData.setFreeR(hkl.getIndex(), (int) data[rFreePlus]);
519 } else if (rFreePlus > 0) {
520 refinementData.setFreeR(hkl.getIndex(), (int) data[rFreePlus]);
521 } else if (rFreeMinus > 0) {
522 refinementData.setFreeR(hkl.getIndex(), (int) data[rFreeMinus]);
523 }
524 }
525 nRead++;
526 } else {
527 HKL tmp = new HKL(ih, ik, il);
528 if (!reflectionList.resolution.inInverseResSqRange(
529 reflectionList.crystal.invressq(tmp))) {
530 nRes++;
531 } else {
532 nIgnore++;
533 }
534 }
535 }
536
537
538 refinementData.generateFsigFfromAnomalousFsigF(anofSigF);
539
540
541 if (logger.isLoggable(Level.INFO)) {
542 sb.append(format(" MTZ file type (machine stamp): %s\n", stampString));
543 sb.append(format(" HKL data is %s\n", transpose ? "transposed" : "not transposed"));
544 sb.append(format(" HKL read in: %d\n", nRead));
545 sb.append(format(" HKL read as friedel mates: %d\n", nFriedel));
546 sb.append(format(" HKL NOT read in (too high resolution): %d\n", nRes));
547 sb.append(format(" HKL NOT read in (not in internal list?): %d\n", nIgnore));
548 sb.append(format(" HKL NOT read in (F/sigF cutoff): %d\n", nCut));
549 sb.append(
550 format(" HKL in internal list: %d", reflectionList.hklList.size()));
551 logger.info(sb.toString());
552 }
553 if (rFree < 0 && rFreePlus < 0 && rFreeMinus < 0) {
554 refinementData.generateRFree();
555 }
556 } catch (EOFException e) {
557 String message = " MTZ end of file reached.";
558 logger.log(Level.WARNING, message, e);
559 return false;
560 } catch (IOException e) {
561 String message = " MTZ IO Exception.";
562 logger.log(Level.WARNING, message, e);
563 return false;
564 }
565
566 return true;
567 }
568
569
570
571
572
573
574
575
576
577
578
579 private boolean readFcs(File mtzFile, ReflectionList reflectionList,
580 DiffractionRefinementData fcData, CompositeConfiguration properties) {
581
582 int nRead, nIgnore, nRes, nFriedel, nCut;
583 ByteOrder byteOrder = ByteOrder.nativeOrder();
584 FileInputStream fileInputStream;
585 DataInputStream dataInputStream;
586
587 StringBuilder sb = new StringBuilder();
588 try {
589 fileInputStream = new FileInputStream(mtzFile);
590 dataInputStream = new DataInputStream(fileInputStream);
591
592 byte[] headerOffset = new byte[4];
593 byte[] bytes = new byte[80];
594 int offset = 0;
595
596
597 dataInputStream.read(bytes, offset, 4);
598
599
600 dataInputStream.read(headerOffset, offset, 4);
601
602
603 dataInputStream.read(bytes, offset, 4);
604 ByteBuffer byteBuffer = ByteBuffer.wrap(bytes);
605 int stamp = byteBuffer.order(ByteOrder.BIG_ENDIAN).getInt();
606 String stampString = Integer.toHexString(stamp);
607 switch (stampString.charAt(0)) {
608 case '1':
609 case '3':
610 if (byteOrder.equals(ByteOrder.LITTLE_ENDIAN)) {
611 byteOrder = ByteOrder.BIG_ENDIAN;
612 }
613 break;
614 case '4':
615 if (byteOrder.equals(ByteOrder.BIG_ENDIAN)) {
616 byteOrder = ByteOrder.LITTLE_ENDIAN;
617 }
618 break;
619 }
620
621 byteBuffer = ByteBuffer.wrap(headerOffset);
622 int headerOffsetI = byteBuffer.order(byteOrder).getInt();
623
624
625 dataInputStream.skipBytes((headerOffsetI - 4) * 4);
626
627 for (Boolean parsing = true; parsing; dataInputStream.read(bytes, offset, 80)) {
628 String mtzString = new String(bytes);
629 parsing = parseHeader(mtzString);
630 }
631
632
633 fc = phiC = fs = phiS = -1;
634 boolean print = true;
635 parseFcColumns(print);
636
637 if (h < 0 || k < 0 || l < 0) {
638 String message = " Fatal error in MTZ file - no H K L indexes?\n";
639 logger.log(Level.SEVERE, message);
640 return false;
641 }
642
643
644 fileInputStream = new FileInputStream(mtzFile);
645 dataInputStream = new DataInputStream(fileInputStream);
646
647
648 dataInputStream.skipBytes(80);
649
650 float[] data = new float[nColumns];
651 HKL mate = new HKL();
652
653
654 ComplexNumber complexNumber = new ComplexNumber();
655 nRead = nIgnore = nRes = nFriedel = nCut = 0;
656 for (int i = 0; i < nReflections; i++) {
657 for (int j = 0; j < nColumns; j++) {
658 dataInputStream.read(bytes, offset, 4);
659 byteBuffer = ByteBuffer.wrap(bytes);
660 data[j] = byteBuffer.order(byteOrder).getFloat();
661 }
662 int ih = (int) data[h];
663 int ik = (int) data[k];
664 int il = (int) data[l];
665 reflectionList.findSymHKL(ih, ik, il, mate, false);
666 HKL hkl = reflectionList.getHKL(mate);
667
668 if (hkl != null) {
669 if (fc > 0 && phiC > 0) {
670 complexNumber.re(data[fc] * cos(toRadians(data[phiC])));
671 complexNumber.im(data[fc] * sin(toRadians(data[phiC])));
672 fcData.setFc(hkl.getIndex(), complexNumber);
673 }
674 if (fs > 0 && phiS > 0) {
675 complexNumber.re(data[fs] * cos(toRadians(data[phiS])));
676 complexNumber.im(data[fs] * sin(toRadians(data[phiS])));
677 fcData.setFs(hkl.getIndex(), complexNumber);
678 }
679 nRead++;
680 } else {
681 HKL tmp = new HKL(ih, ik, il);
682 if (!reflectionList.resolution.inInverseResSqRange(
683 reflectionList.crystal.invressq(tmp))) {
684 nRes++;
685 } else {
686 nIgnore++;
687 }
688 }
689 }
690
691 if (logger.isLoggable(Level.INFO)) {
692 sb.append(format(" MTZ file type (machine stamp): %s\n", stampString));
693 sb.append(format(" Fc HKL read in: %d\n", nRead));
694 sb.append(format(" Fc HKL read as friedel mates: %d\n", nFriedel));
695 sb.append(format(" Fc HKL NOT read in (too high resolution): %d\n", nRes));
696 sb.append(format(" Fc HKL NOT read in (not in internal list?): %d\n", nIgnore));
697 sb.append(format(" Fc HKL NOT read in (F/sigF cutoff): %d\n", nCut));
698 sb.append(
699 format(
700 " HKL in internal list: %d\n",
701 reflectionList.hklList.size()));
702 logger.info(sb.toString());
703 }
704 } catch (EOFException e) {
705 String message = " MTZ end of file reached.";
706 logger.log(Level.WARNING, message, e);
707 return false;
708 } catch (IOException e) {
709 String message = " MTZ IO Exception.";
710 logger.log(Level.WARNING, message, e);
711 return false;
712 }
713
714 return true;
715 }
716
717
718
719
720
721
722
723 private Boolean parseHeader(String str) {
724 boolean parsing = true;
725 Dataset dataSet;
726 String[] strArray = str.split("\\s+");
727
728 if (headerParsed) {
729 return Header.toHeader(strArray[0]) != Header.END;
730 }
731
732 switch (Header.toHeader(strArray[0])) {
733 case TITLE:
734 title = str.substring(5);
735 break;
736 case NCOL:
737 nColumns = parseInt(strArray[1]);
738 nReflections = parseInt(strArray[2]);
739 int nBatches = parseInt(strArray[3]);
740 break;
741 case SYMINF:
742 String[] tmp = str.split("\'+");
743 spaceGroupNum = parseInt(strArray[4]);
744 if (tmp.length > 1) {
745 spaceGroupName = tmp[1];
746 }
747 break;
748 case RESO:
749 double r1 = sqrt(1.0 / parseFloat(strArray[1]));
750 double r2 = sqrt(1.0 / parseFloat(strArray[2]));
751 resLow = max(r1, r2);
752 resHigh = min(r1, r2);
753 break;
754 case NDIF:
755 int ndif = parseInt(strArray[1]);
756 break;
757 case COL:
758 case COLUMN:
759 int nDataSet = parseInt(strArray[5]);
760 if (nDataSet == 0) {
761 dsetOffset = 0;
762 }
763 Column column = new Column();
764 columns.add(column);
765 column.label = strArray[1];
766 column.type = strArray[2].charAt(0);
767 column.id = nDataSet;
768 column.min = parseDouble(strArray[3]);
769 column.max = parseDouble(strArray[4]);
770 break;
771 case PROJECT:
772 nDataSet = parseInt(strArray[1]);
773 if (nDataSet == 0) {
774 dsetOffset = 0;
775 }
776 try {
777 dataSet = dataSets.get(nDataSet - dsetOffset);
778 } catch (IndexOutOfBoundsException e) {
779 dataSet = new Dataset();
780 dataSets.add(dataSet);
781 }
782 dataSet.project = strArray[2];
783 break;
784 case DATASET:
785 nDataSet = parseInt(strArray[1]);
786 if (nDataSet == 0) {
787 dsetOffset = 0;
788 }
789 try {
790 dataSet = dataSets.get(nDataSet - dsetOffset);
791 } catch (IndexOutOfBoundsException e) {
792 dataSet = new Dataset();
793 dataSets.add(dataSet);
794 }
795 dataSet.dataset = strArray[2];
796 break;
797 case DCELL:
798 nDataSet = Integer.parseInt(strArray[1]);
799 if (nDataSet == 0) {
800 dsetOffset = 0;
801 }
802 try {
803 dataSet = dataSets.get(nDataSet - dsetOffset);
804 } catch (IndexOutOfBoundsException e) {
805 dataSet = new Dataset();
806 dataSets.add(dataSet);
807 }
808 dataSet.cell[0] = parseDouble(strArray[2]);
809 dataSet.cell[1] = parseDouble(strArray[3]);
810 dataSet.cell[2] = parseDouble(strArray[4]);
811 dataSet.cell[3] = parseDouble(strArray[5]);
812 dataSet.cell[4] = parseDouble(strArray[6]);
813 dataSet.cell[5] = parseDouble(strArray[7]);
814 break;
815 case DWAVEL:
816 nDataSet = parseInt(strArray[1]);
817 if (nDataSet == 0) {
818 dsetOffset = 0;
819 }
820 try {
821 dataSet = dataSets.get(nDataSet - dsetOffset);
822 } catch (IndexOutOfBoundsException e) {
823 dataSet = new Dataset();
824 dataSets.add(dataSet);
825 }
826 dataSet.lambda = parseDouble(strArray[2]);
827 break;
828 case END:
829 headerParsed = true;
830 parsing = false;
831 break;
832 default:
833 break;
834 }
835
836 return parsing;
837 }
838
839
840
841
842
843
844 private void parseColumns(boolean print) {
845
846 int nc = 0;
847 StringBuilder sb = new StringBuilder();
848 for (Iterator<Column> i = columns.iterator(); i.hasNext(); nc++) {
849 Column column = i.next();
850 String label = column.label.trim();
851 if (label.equalsIgnoreCase("H") && column.type == 'H') {
852 h = nc;
853 } else if (label.equalsIgnoreCase("K") && column.type == 'H') {
854 k = nc;
855 } else if (label.equalsIgnoreCase("L") && column.type == 'H') {
856 l = nc;
857 } else if ((label.equalsIgnoreCase("free")
858 || label.equalsIgnoreCase("freer")
859 || label.equalsIgnoreCase("freerflag")
860 || label.equalsIgnoreCase("freer_flag")
861 || label.equalsIgnoreCase("rfree")
862 || label.equalsIgnoreCase("rfreeflag")
863 || label.equalsIgnoreCase("r-free-flags")
864 || label.equalsIgnoreCase("test")
865 || Strings.CI.equals(label, rFreeString))
866 && column.type == 'I') {
867 sb.append(format(" Reading R Free column: \"%s\"\n", column.label));
868 rFree = nc;
869 } else if ((label.equalsIgnoreCase("free(+)")
870 || label.equalsIgnoreCase("freer(+)")
871 || label.equalsIgnoreCase("freerflag(+)")
872 || label.equalsIgnoreCase("freer_flag(+)")
873 || label.equalsIgnoreCase("rfree(+)")
874 || label.equalsIgnoreCase("rfreeflag(+)")
875 || label.equalsIgnoreCase("r-free-flags(+)")
876 || label.equalsIgnoreCase("test(+)")
877 || Strings.CI.equals(label + "(+)", rFreeString))
878 && column.type == 'I') {
879 rFreePlus = nc;
880 } else if ((label.equalsIgnoreCase("free(-)")
881 || label.equalsIgnoreCase("freer(-)")
882 || label.equalsIgnoreCase("freerflag(-)")
883 || label.equalsIgnoreCase("freer_flag(-)")
884 || label.equalsIgnoreCase("rfree(-)")
885 || label.equalsIgnoreCase("rfreeflag(-)")
886 || label.equalsIgnoreCase("r-free-flags(-)")
887 || label.equalsIgnoreCase("test(-)")
888 || Strings.CI.equals(label + "(-)", rFreeString))
889 && column.type == 'I') {
890 rFreeMinus = nc;
891 } else if ((label.equalsIgnoreCase("f")
892 || label.equalsIgnoreCase("fp")
893 || label.equalsIgnoreCase("fo")
894 || label.equalsIgnoreCase("fobs")
895 || label.equalsIgnoreCase("f-obs")
896 || Strings.CI.equals(label, foString))
897 && column.type == 'F') {
898 sb.append(format(" Reading Fo column: \"%s\"\n", column.label));
899 fo = nc;
900 } else if ((label.equalsIgnoreCase("f(+)")
901 || label.equalsIgnoreCase("fp(+)")
902 || label.equalsIgnoreCase("fo(+)")
903 || label.equalsIgnoreCase("fobs(+)")
904 || label.equalsIgnoreCase("f-obs(+)")
905 || Strings.CI.equals(label + "(+)", foString))
906 && column.type == 'G') {
907 fPlus = nc;
908 } else if ((label.equalsIgnoreCase("f(-)")
909 || label.equalsIgnoreCase("fp(-)")
910 || label.equalsIgnoreCase("fo(-)")
911 || label.equalsIgnoreCase("fobs(-)")
912 || label.equalsIgnoreCase("f-obs(-)")
913 || Strings.CI.equals(label + "(-)", foString))
914 && column.type == 'G') {
915 fMinus = nc;
916 } else if ((label.equalsIgnoreCase("sigf")
917 || label.equalsIgnoreCase("sigfp")
918 || label.equalsIgnoreCase("sigfo")
919 || label.equalsIgnoreCase("sigfobs")
920 || label.equalsIgnoreCase("sigf-obs")
921 || Strings.CI.equals(label, sigFoString))
922 && column.type == 'Q') {
923 sb.append(format(" Reading sigFo column: \"%s\"\n", column.label));
924 sigFo = nc;
925 } else if ((label.equalsIgnoreCase("sigf(+)")
926 || label.equalsIgnoreCase("sigfp(+)")
927 || label.equalsIgnoreCase("sigfo(+)")
928 || label.equalsIgnoreCase("sigfobs(+)")
929 || label.equalsIgnoreCase("sigf-obs(+)")
930 || Strings.CI.equals(label + "(+)", sigFoString))
931 && column.type == 'L') {
932 sigFPlus = nc;
933 } else if ((label.equalsIgnoreCase("sigf(-)")
934 || label.equalsIgnoreCase("sigfp(-)")
935 || label.equalsIgnoreCase("sigfo(-)")
936 || label.equalsIgnoreCase("sigfobs(-)")
937 || label.equalsIgnoreCase("sigf-obs(-)")
938 || Strings.CI.equals(label + "(-)", sigFoString))
939 && column.type == 'L') {
940 sigFMinus = nc;
941 }
942 }
943 if (fo < 0 && sigFo < 0 && fPlus > 0 && sigFPlus > 0 && fMinus > 0 && sigFMinus > 0) {
944 sb.append(" Reading Fplus/Fminus column to fill in Fo\n");
945 }
946 if (logger.isLoggable(Level.INFO) && print) {
947 logger.info(sb.toString());
948 }
949 }
950
951
952
953
954
955
956 private void parseFcColumns(boolean print) {
957 int nc = 0;
958 StringBuilder sb = new StringBuilder();
959 for (Iterator<Column> i = columns.iterator(); i.hasNext(); nc++) {
960 Column column = i.next();
961 String label = column.label.trim();
962 if (label.equalsIgnoreCase("H") && column.type == 'H') {
963 h = nc;
964 } else if (label.equalsIgnoreCase("K") && column.type == 'H') {
965 k = nc;
966 } else if (label.equalsIgnoreCase("L") && column.type == 'H') {
967 l = nc;
968 } else if ((label.equalsIgnoreCase("fc") || label.equalsIgnoreCase("fcalc"))
969 && column.type == 'F') {
970 sb.append(format(" Reading Fc column: \"%s\"\n", column.label));
971 fc = nc;
972 } else if ((label.equalsIgnoreCase("phic")
973 || label.equalsIgnoreCase("phifc")
974 || label.equalsIgnoreCase("phicalc")
975 || label.equalsIgnoreCase("phifcalc"))
976 && column.type == 'P') {
977 sb.append(format(" Reading phiFc column: \"%s\"\n", column.label));
978 phiC = nc;
979 } else if ((label.equalsIgnoreCase("fs") || label.equalsIgnoreCase("fscalc"))
980 && column.type == 'F') {
981 sb.append(format(" Reading Fs column: \"%s\"\n", column.label));
982 fs = nc;
983 } else if ((label.equalsIgnoreCase("phis")
984 || label.equalsIgnoreCase("phifs")
985 || label.equalsIgnoreCase("phiscalc")
986 || label.equalsIgnoreCase("phifscalc"))
987 && column.type == 'P') {
988 sb.append(format(" Reading phiFs column: \"%s\"\n", column.label));
989 phiS = nc;
990 }
991 }
992 if (logger.isLoggable(Level.INFO) && print) {
993 logger.info(sb.toString());
994 }
995 }
996
997
998
999
1000
1001
1002 int getnColumns() {
1003 return nColumns;
1004 }
1005
1006
1007
1008
1009
1010
1011 int getnReflections() {
1012 return nReflections;
1013 }
1014
1015
1016
1017
1018
1019
1020 int getSpaceGroupNum() {
1021 return spaceGroupNum;
1022 }
1023
1024 private enum Header {
1025 VERS,
1026 TITLE,
1027 NCOL,
1028 SORT,
1029 SYMINF,
1030 SYMM,
1031 RESO,
1032 VALM,
1033 COL,
1034 COLUMN,
1035 NDIF,
1036 PROJECT,
1037 CRYSTAL,
1038 DATASET,
1039 DCELL,
1040 DWAVEL,
1041 BATCH,
1042 END,
1043 NOVALUE;
1044
1045 public static Header toHeader(String str) {
1046 try {
1047 return valueOf(str);
1048 } catch (Exception ex) {
1049 return NOVALUE;
1050 }
1051 }
1052 }
1053
1054 private static class Column {
1055
1056 public String label;
1057 public char type;
1058 public int id;
1059 public double min, max;
1060 }
1061
1062 private static class Dataset {
1063
1064 public double lambda;
1065 public double[] cell = new double[6];
1066 String project;
1067 String dataset;
1068 }
1069 }