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