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