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