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 ffx.utilities.TinkerUtils.version;
41 import static java.lang.Double.isNaN;
42 import static java.lang.String.format;
43 import static org.apache.commons.math3.util.FastMath.max;
44 import static org.apache.commons.math3.util.FastMath.min;
45 import static org.apache.commons.math3.util.FastMath.toDegrees;
46
47 import ffx.crystal.Crystal;
48 import ffx.crystal.HKL;
49 import ffx.crystal.ReflectionList;
50 import ffx.crystal.ReflectionSpline;
51 import ffx.crystal.SpaceGroup;
52 import ffx.crystal.SymOp;
53 import ffx.xray.DiffractionRefinementData;
54 import java.io.DataOutputStream;
55 import java.io.File;
56 import java.io.FileOutputStream;
57 import java.nio.ByteBuffer;
58 import java.nio.ByteOrder;
59 import java.text.SimpleDateFormat;
60 import java.util.Date;
61 import java.util.Vector;
62 import java.util.logging.Level;
63 import java.util.logging.Logger;
64
65
66
67
68
69
70
71
72
73
74 public class MTZWriter {
75
76 private static final Logger logger = Logger.getLogger(MTZWriter.class.getName());
77 private final String fileName;
78 private final ReflectionList reflectionList;
79 private final DiffractionRefinementData refinementData;
80 private final Crystal crystal;
81 private final SpaceGroup spaceGroup;
82 private final ReflectionSpline spline;
83 private final int n;
84 private final int nCol;
85 private final int mtzType;
86
87
88
89
90
91
92
93
94 public MTZWriter(
95 ReflectionList reflectionList, DiffractionRefinementData refinementData, String filename) {
96 this(reflectionList, refinementData, filename, MTZType.ALL);
97 }
98
99
100
101
102
103
104
105
106
107 public MTZWriter(
108 ReflectionList reflectionList,
109 DiffractionRefinementData refinementData,
110 String filename,
111 int mtzType) {
112 this.reflectionList = reflectionList;
113 this.refinementData = refinementData;
114 this.crystal = reflectionList.crystal;
115 this.spaceGroup = crystal.spaceGroup;
116 this.spline = new ReflectionSpline(reflectionList, refinementData.spline.length);
117 this.fileName = filename;
118
119
120 this.n = refinementData.n - 1;
121 this.mtzType = mtzType;
122 switch (mtzType) {
123 case MTZType.DATAONLY:
124 this.nCol = 6;
125 break;
126 case MTZType.FCONLY:
127 this.nCol = 7;
128 break;
129 default:
130 this.nCol = 18;
131 break;
132 }
133 }
134
135
136 public void write() {
137 ByteOrder byteOrder = ByteOrder.nativeOrder();
138 FileOutputStream fileOutputStream;
139 DataOutputStream dataOutputStream;
140
141 File file = version(new File(fileName));
142 String name = file.getName();
143
144 try {
145 if (logger.isLoggable(Level.INFO)) {
146 logger.info(format("\n Saving %s", name));
147 }
148
149 fileOutputStream = new FileOutputStream(file);
150 dataOutputStream = new DataOutputStream(fileOutputStream);
151
152 byte[] bytes = new byte[80];
153 int offset = 0;
154 int writeLen = 0;
155
156 int iMapData;
157 float fMapData;
158
159
160 StringBuilder sb = new StringBuilder();
161 sb.append("MTZ ");
162 dataOutputStream.writeBytes(sb.toString());
163
164
165 int headerOffset = n * nCol + 21;
166 ByteBuffer byteBuffer = ByteBuffer.wrap(bytes);
167 byteBuffer.order(byteOrder).putInt(headerOffset);
168
169
170
171 if (ByteOrder.nativeOrder().equals(ByteOrder.LITTLE_ENDIAN)) {
172 iMapData = 0x4441;
173 } else {
174 iMapData = 0x1111;
175 }
176 byteBuffer.order(byteOrder).putInt(iMapData);
177 dataOutputStream.write(bytes, offset, 8);
178
179 sb = new StringBuilder();
180 sb.append(" ");
181 sb.setLength(68);
182 dataOutputStream.writeBytes(sb.toString());
183
184
185 Vector<String> colname = new Vector<>(nCol);
186 char[] colType = new char[nCol];
187 double[] res = new double[2];
188 res[0] = Double.POSITIVE_INFINITY;
189 res[1] = Double.NEGATIVE_INFINITY;
190 float[][] colMinMax = new float[nCol][2];
191 for (int i = 0; i < nCol; i++) {
192 colMinMax[i][0] = Float.POSITIVE_INFINITY;
193 colMinMax[i][1] = Float.NEGATIVE_INFINITY;
194 }
195 ReflectionSpline sigmaASpline =
196 new ReflectionSpline(reflectionList, refinementData.sigmaA.length);
197 int col = 0;
198
199 colname.add("H");
200 colType[col++] = 'H';
201 colname.add("K");
202 colType[col++] = 'H';
203 colname.add("L");
204 colType[col++] = 'H';
205 writeLen += 12;
206
207 if (mtzType != MTZType.FCONLY) {
208 colname.add("FO");
209 colType[col++] = 'F';
210 colname.add("SIGFO");
211 colType[col++] = 'Q';
212 colname.add("FreeR");
213 colType[col++] = 'I';
214 writeLen += 12;
215 }
216
217 if (mtzType != MTZType.DATAONLY) {
218 colname.add("Fs");
219 colType[col++] = 'F';
220 colname.add("PHIFs");
221 colType[col++] = 'P';
222 colname.add("Fc");
223 colType[col++] = 'F';
224 colname.add("PHIFc");
225 colType[col++] = 'P';
226 writeLen += 16;
227 }
228
229 if (mtzType == MTZType.ALL) {
230 colname.add("FOM");
231 colType[col++] = 'W';
232 colname.add("PHIW");
233 colType[col++] = 'P';
234 colname.add("SigmaAs");
235 colType[col++] = 'F';
236 colname.add("SigmaAw");
237 colType[col++] = 'Q';
238 colname.add("FWT");
239 colType[col++] = 'F';
240 colname.add("PHWT");
241 colType[col++] = 'P';
242 colname.add("DELFWT");
243 colType[col++] = 'F';
244 colname.add("PHDELWT");
245 colType[col] = 'P';
246 writeLen += 32;
247 }
248
249 for (HKL ih : reflectionList.hklList) {
250 col = 0;
251 int i = ih.getIndex();
252
253
254 if (ih.getH() == 0 && ih.getK() == 0 && ih.getL() == 0) {
255 continue;
256 }
257
258 double ss = crystal.invressq(ih);
259 res[0] = min(ss, res[0]);
260 res[1] = max(ss, res[1]);
261
262
263 fMapData = ih.getH();
264 colMinMax[col][0] = min(fMapData, colMinMax[0][0]);
265 colMinMax[col][1] = max(fMapData, colMinMax[0][1]);
266 byteBuffer.rewind();
267 byteBuffer.order(byteOrder).putFloat(fMapData);
268 col++;
269
270 fMapData = ih.getK();
271 colMinMax[col][0] = min(fMapData, colMinMax[1][0]);
272 colMinMax[col][1] = max(fMapData, colMinMax[1][1]);
273 byteBuffer.order(byteOrder).putFloat(fMapData);
274 col++;
275
276 fMapData = ih.getL();
277 colMinMax[col][0] = min(fMapData, colMinMax[2][0]);
278 colMinMax[col][1] = max(fMapData, colMinMax[2][1]);
279 byteBuffer.order(byteOrder).putFloat(fMapData);
280 col++;
281
282 if (mtzType != MTZType.FCONLY) {
283
284 fMapData = (float) refinementData.getF(i);
285 if (!isNaN(fMapData)) {
286 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
287 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
288 }
289 byteBuffer.order(byteOrder).putFloat(fMapData);
290 col++;
291 fMapData = (float) refinementData.getSigF(i);
292 if (!isNaN(fMapData)) {
293 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
294 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
295 }
296 byteBuffer.order(byteOrder).putFloat(fMapData);
297 col++;
298
299
300 fMapData = (float) refinementData.getFreeR(i);
301 if (!isNaN(fMapData)) {
302 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
303 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
304 }
305 byteBuffer.order(byteOrder).putFloat(fMapData);
306 col++;
307 }
308
309 if (mtzType == MTZType.FCONLY) {
310
311 fMapData = (float) refinementData.fsF(i);
312 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
313 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
314 byteBuffer.order(byteOrder).putFloat(fMapData);
315 col++;
316 fMapData = (float) toDegrees(refinementData.fsPhi(i));
317 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
318 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
319 byteBuffer.order(byteOrder).putFloat(fMapData);
320 col++;
321
322
323 fMapData = (float) refinementData.fcF(i);
324 if (!isNaN(fMapData)) {
325 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
326 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
327 }
328 byteBuffer.order(byteOrder).putFloat(fMapData);
329 col++;
330 fMapData = (float) Math.toDegrees(refinementData.fcPhi(i));
331 if (!isNaN(fMapData)) {
332 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
333 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
334 }
335 byteBuffer.order(byteOrder).putFloat(fMapData);
336 col++;
337 }
338
339 if (mtzType == MTZType.ALL) {
340
341 fMapData = (float) refinementData.fsF(i);
342 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
343 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
344 byteBuffer.order(byteOrder).putFloat(fMapData);
345 col++;
346 fMapData = (float) toDegrees(refinementData.fsPhi(i));
347 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
348 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
349 byteBuffer.order(byteOrder).putFloat(fMapData);
350 col++;
351
352
353 fMapData = (float) refinementData.fcTotF(i);
354 if (!isNaN(fMapData)) {
355 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
356 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
357 }
358 byteBuffer.order(byteOrder).putFloat(fMapData);
359 col++;
360 fMapData = (float) toDegrees(refinementData.fcTotPhi(i));
361 if (!isNaN(fMapData)) {
362 colMinMax[col][0] = Math.min(fMapData, colMinMax[col][0]);
363 colMinMax[col][1] = Math.max(fMapData, colMinMax[col][1]);
364 }
365 byteBuffer.order(byteOrder).putFloat(fMapData);
366 col++;
367
368
369 fMapData = (float) refinementData.fomPhi[i][0];
370 if (!isNaN(fMapData)) {
371 colMinMax[col][0] = Math.min(fMapData, colMinMax[col][0]);
372 colMinMax[col][1] = Math.max(fMapData, colMinMax[col][1]);
373 }
374 byteBuffer.order(byteOrder).putFloat(fMapData);
375 col++;
376 fMapData = (float) toDegrees(refinementData.fomPhi[i][1]);
377 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
378 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
379 byteBuffer.order(byteOrder).putFloat(fMapData);
380 col++;
381
382
383 spline.f(ss, refinementData.spline);
384 double sa = sigmaASpline.f(ss, refinementData.sigmaA);
385 double wa = sigmaASpline.f(ss, refinementData.sigmaW);
386
387
388 fMapData = (float) sa;
389 if (!isNaN(fMapData)) {
390 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
391 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
392 }
393 byteBuffer.order(byteOrder).putFloat(fMapData);
394 col++;
395 fMapData = (float) wa;
396 if (!isNaN(fMapData)) {
397 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
398 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
399 }
400 byteBuffer.order(byteOrder).putFloat(fMapData);
401 col++;
402
403
404 fMapData = (float) refinementData.FoFc2F(i);
405 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
406 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
407 byteBuffer.order(byteOrder).putFloat(fMapData);
408 col++;
409 fMapData = (float) toDegrees(refinementData.FoFc2Phi(i));
410 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
411 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
412 byteBuffer.order(byteOrder).putFloat(fMapData);
413 col++;
414 fMapData = (float) refinementData.foFc1F(i);
415 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
416 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
417 byteBuffer.order(byteOrder).putFloat(fMapData);
418 col++;
419 fMapData = (float) toDegrees(refinementData.foFc1Phi(i));
420 colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
421 colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
422 byteBuffer.order(byteOrder).putFloat(fMapData);
423 }
424
425 dataOutputStream.write(bytes, offset, writeLen);
426 }
427
428
429 sb = new StringBuilder();
430 sb.append("VERS MTZ:V1.1 ");
431 while (sb.length() < 80) {
432 sb.append(" ");
433 }
434 dataOutputStream.writeBytes(sb.toString());
435
436 Date now = new Date();
437 SimpleDateFormat sdf = new SimpleDateFormat("yyyy-MM-dd'T'HH:mm:ss ");
438 sb = new StringBuilder();
439 sb.append("TITLE FFX output: ").append(sdf.format(now));
440 while (sb.length() < 80) {
441 sb.append(" ");
442 }
443 dataOutputStream.writeBytes(sb.toString());
444
445 sb = new StringBuilder();
446 sb.append(String.format("NCOL %8d %12d %8d", nCol, n, 0));
447 while (sb.length() < 80) {
448 sb.append(" ");
449 }
450 dataOutputStream.writeBytes(sb.toString());
451
452 sb = new StringBuilder();
453 sb.append("SORT 0 0 0 0 0 ");
454 while (sb.length() < 80) {
455 sb.append(" ");
456 }
457 dataOutputStream.writeBytes(sb.toString());
458
459 sb = new StringBuilder();
460 char cdata = spaceGroup.shortName.charAt(0);
461 if (cdata == 'H') {
462 cdata = 'R';
463 }
464 sb.append(
465 String.format(
466 "SYMINF %3d %2d %c %5d %22s %5s",
467 spaceGroup.getNumberOfSymOps(),
468 spaceGroup.numPrimitiveSymEquiv,
469 cdata,
470 spaceGroup.number,
471 "'" + spaceGroup.shortName + "'",
472 spaceGroup.pointGroupName));
473 while (sb.length() < 80) {
474 sb.append(" ");
475 }
476 dataOutputStream.writeBytes(sb.toString());
477
478 for (int i = 0; i < spaceGroup.symOps.size(); i++) {
479 sb = new StringBuilder();
480 sb.append("SYMM ");
481 SymOp symop = spaceGroup.symOps.get(i);
482 sb.append(symop.toXYZString());
483 while (sb.length() < 80) {
484 sb.append(" ");
485 }
486 dataOutputStream.writeBytes(sb.toString());
487 }
488
489 sb = new StringBuilder();
490 sb.append(String.format("RESO %8.6f%13s%8.6f", res[0], " ", res[1]));
491 while (sb.length() < 80) {
492 sb.append(" ");
493 }
494 dataOutputStream.writeBytes(sb.toString());
495
496 sb = new StringBuilder();
497 sb.append("VALM NAN ");
498 while (sb.length() < 80) {
499 sb.append(" ");
500 }
501 dataOutputStream.writeBytes(sb.toString());
502
503 sb = new StringBuilder();
504 sb.append("NDIF 1 ");
505 while (sb.length() < 80) {
506 sb.append(" ");
507 }
508 dataOutputStream.writeBytes(sb.toString());
509
510 sb = new StringBuilder();
511 sb.append("PROJECT 1 project ");
512 while (sb.length() < 80) {
513 sb.append(" ");
514 }
515 dataOutputStream.writeBytes(sb.toString());
516
517 sb = new StringBuilder();
518 sb.append("CRYSTAL 1 crystal ");
519 while (sb.length() < 80) {
520 sb.append(" ");
521 }
522 dataOutputStream.writeBytes(sb.toString());
523
524 sb = new StringBuilder();
525 sb.append("DATASET 1 dataset ");
526 while (sb.length() < 80) {
527 sb.append(" ");
528 }
529 dataOutputStream.writeBytes(sb.toString());
530
531 for (int j = 0; j < nCol; j++) {
532 sb = new StringBuilder();
533 sb.append(
534 String.format(
535 "COLUMN %-30s %c %17.4f %17.4f 1",
536 colname.get(j), colType[j], colMinMax[j][0], colMinMax[j][1]));
537 dataOutputStream.writeBytes(sb.toString());
538 }
539
540 sb = new StringBuilder();
541 sb.append(
542 String.format(
543 "CELL %10.4f %9.4f %9.4f %9.4f %9.4f %9.4f ",
544 crystal.a, crystal.b, crystal.c, crystal.alpha, crystal.beta, crystal.gamma));
545 while (sb.length() < 80) {
546 sb.append(" ");
547 }
548 dataOutputStream.writeBytes(sb.toString());
549
550 sb = new StringBuilder();
551 sb.append(
552 String.format(
553 "DCELL %9d %10.4f %9.4f %9.4f %9.4f %9.4f %9.4f ",
554 1, crystal.a, crystal.b, crystal.c, crystal.alpha, crystal.beta, crystal.gamma));
555 while (sb.length() < 80) {
556 sb.append(" ");
557 }
558 dataOutputStream.writeBytes(sb.toString());
559
560 sb = new StringBuilder();
561 sb.append("DWAVEL 1 1.00000 ");
562 while (sb.length() < 80) {
563 sb.append(" ");
564 }
565 dataOutputStream.writeBytes(sb.toString());
566
567 sb = new StringBuilder();
568 sb.append("END ");
569 while (sb.length() < 80) {
570 sb.append(" ");
571 }
572 dataOutputStream.writeBytes(sb.toString());
573
574 sb = new StringBuilder();
575 sb.append("MTZENDOFHEADERS ");
576 while (sb.length() < 80) {
577 sb.append(" ");
578 }
579 dataOutputStream.writeBytes(sb.toString());
580 dataOutputStream.close();
581
582 if (logger.isLoggable(Level.INFO)) {
583 logger.info(format(" Wrote MTZ to file %s", file.getAbsolutePath()));
584 }
585
586 } catch (Exception e) {
587 String message = "Fatal exception evaluating structure factors.\n";
588 logger.log(Level.SEVERE, message, e);
589 }
590 }
591
592
593 public interface MTZType {
594
595
596 int DATAONLY = 1;
597
598 int FCONLY = 2;
599
600 int ALL = 3;
601 }
602 }