View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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   * This class parses CCP4 MTZ files.<br>
76   *
77   * @author Timothy D. Fenn<br>
78   * @see <a href="http://www.ccp4.ac.uk/html/maplib.html" target="_blank">CCP4 map format</a>
79   * @see <a href="http://www.ccp4.ac.uk/dist/html/library.html" target="_blank">CCP4 library
80   *     documentation</a>
81   * @since 1.0
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   /** Constructor for MTZFilter. */
103   public MTZFilter() {
104   }
105 
106   /**
107    * Average the computed structure factors for two systems.
108    *
109    * @param mtzFile1 This file will be overwritten and become the new average.
110    * @param mtzFile2 Second MTZ file.
111    * @param reflectionlist List of HKLs.
112    * @param iter The iteration in the running average.
113    * @param properties The CompositeConfiguration defines the properties of each system.
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     // compute running average using mtzFile1 as current average
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     // overwrite original MTZ
139     MTZWriter mtzOut = new MTZWriter(reflectionlist, fcdata1, mtzFile1.getName(), MTZType.FCONLY);
140     mtzOut.write();
141   }
142 
143   /** {@inheritDoc} */
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       // Eat "MTZ" title.
158       dataInputStream.read(bytes, offset, 4);
159 
160       // Header offset.
161       dataInputStream.read(headerOffset, offset, 4);
162 
163       // Machine stamp.
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       // skip to header and parse
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     // column identifiers
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   /** {@inheritDoc} */
275   @Override
276   public ReflectionList getReflectionList(File mtzFile) {
277     return getReflectionList(mtzFile, null);
278   }
279 
280   /** {@inheritDoc} */
281   @Override
282   public double getResolution(File mtzFile, Crystal crystal) {
283     ReflectionList reflectionList = getReflectionList(mtzFile, null);
284     return reflectionList.getMaxResolution();
285   }
286 
287   /** printHeader */
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   /** {@inheritDoc} */
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       // Eat "MTZ" title.
349       dataInputStream.read(bytes, offset, 4);
350       // String mtzstr;
351 
352       // Header offset.
353       dataInputStream.read(headerOffset, offset, 4);
354 
355       // Machine stamp.
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       // skip to header and parse
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       // column identifiers
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       // Reopen to start at beginning.
404       fileInputStream = new FileInputStream(mtzFile);
405       dataInputStream = new DataInputStream(fileInputStream);
406 
407       // Skip initial header.
408       dataInputStream.skipBytes(80);
409 
410       // Check if HKLs need to be transposed or not.
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       // Reopen to start at beginning
479       fileInputStream = new FileInputStream(mtzFile);
480       dataInputStream = new DataInputStream(fileInputStream);
481 
482       // Skip initial header
483       dataInputStream.skipBytes(80);
484 
485       // Read in data
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               // not sure what the correct thing to do here is?
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       // Set up fsigf from F+ and F-.
565       refinementData.generateFsigFfromAnomalousFsigF(anofSigF);
566 
567       // Log results.
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    * Read the structure factors.
598    *
599    * @param mtzFile a {@link java.io.File} object.
600    * @param reflectionList a {@link ffx.crystal.ReflectionList} object.
601    * @param fcData a {@link ffx.xray.DiffractionRefinementData} object.
602    * @param properties a {@link org.apache.commons.configuration2.CompositeConfiguration}
603    *     object.
604    * @return a boolean.
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       // Eat "MTZ" title.
627       dataInputStream.read(bytes, offset, 4);
628 
629       // Header offset.
630       dataInputStream.read(headerOffset, offset, 4);
631 
632       // Machine stamp.
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       // Skip to header and parse.
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       // Column identifiers.
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       // Reopen to start at beginning.
674       fileInputStream = new FileInputStream(mtzFile);
675       dataInputStream = new DataInputStream(fileInputStream);
676 
677       // Skip initial header.
678       dataInputStream.skipBytes(80);
679 
680       float[] data = new float[nColumns];
681       HKL mate = new HKL();
682 
683       // Read in data.
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    * Parse the header.
749    *
750    * @param str The header to parse.
751    * @return
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    * Parse columns.
871    *
872    * @param print
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    * Parse Fc columns.
983    *
984    * @param print
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    * Getter for the field <code>nColumns</code>.
1029    *
1030    * @return the nColumns
1031    */
1032   int getnColumns() {
1033     return nColumns;
1034   }
1035 
1036   /**
1037    * Getter for the field <code>nReflections</code>.
1038    *
1039    * @return the nReflections
1040    */
1041   int getnReflections() {
1042     return nReflections;
1043   }
1044 
1045   /**
1046    * Getter for the field <code>spaceGroupNum</code>.
1047    *
1048    * @return the spaceGroupNum
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 }