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-2023.
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.crystal;
39  
40  import static ffx.utilities.TinkerUtils.version;
41  import static java.lang.String.format;
42  import static org.apache.commons.math3.util.FastMath.pow;
43  import static org.apache.commons.math3.util.FastMath.sqrt;
44  
45  import java.io.DataOutputStream;
46  import java.io.File;
47  import java.io.FileOutputStream;
48  import java.nio.ByteBuffer;
49  import java.nio.ByteOrder;
50  import java.util.logging.Level;
51  import java.util.logging.Logger;
52  
53  /**
54   * CCP4MapWriter class.
55   *
56   * @author Timothy D. Fenn
57   * @see <ul>
58   *     <li><a href="http://www.ccp4.ac.uk/html/maplib.html" target="_blank">CCP4 map format</a>
59   *     <li><a href="http://www.ccp4.ac.uk/dist/html/library.html" target="_blank">CCP4 library
60   *     documentation </a>
61   *     <li><a href="http://www.ccp4.ac.uk/html/maplib.html" target="_blank">CCP4 map format </a>
62   *     <li><a href="http://www.ccp4.ac.uk/dist/html/library.html" target="_blank">CCP4 library
63   *     documentation </a>
64   *     </ul>
65   * @since 1.0
66   */
67  public class CCP4MapWriter {
68  
69    private static final Logger logger = Logger.getLogger(CCP4MapWriter.class.getName());
70    private final String filename;
71    private final Crystal crystal;
72    private final int extx, exty, extz;
73    private final int orix, oriy, oriz, nx, ny, nz;
74    private int stride;
75  
76    /**
77     * construct mapwriter object
78     *
79     * @param extx slices in x
80     * @param exty slices in y
81     * @param extz slices in z
82     * @param crystal {@link Crystal} object
83     * @param filename output filename
84     */
85    public CCP4MapWriter(int extx, int exty, int extz, Crystal crystal, String filename) {
86      this.orix = 0;
87      this.oriy = 0;
88      this.oriz = 0;
89      this.extx = extx;
90      this.exty = exty;
91      this.extz = extz;
92      this.nx = extx;
93      this.ny = exty;
94      this.nz = extz;
95      this.crystal = crystal;
96      this.filename = filename;
97      this.stride = 2;
98    }
99  
100   /**
101    * Constructor for CCP4MapWriter.
102    *
103    * @param orix an int.
104    * @param oriy an int.
105    * @param oriz an int.
106    * @param extx an int.
107    * @param exty an int.
108    * @param extz an int.
109    * @param nx an int.
110    * @param ny an int.
111    * @param nz an int.
112    * @param crystal a {@link ffx.crystal.Crystal} object.
113    * @param filename a {@link java.lang.String} object.
114    */
115   public CCP4MapWriter(
116       int orix, int oriy, int oriz,
117       int extx, int exty, int extz,
118       int nx, int ny, int nz,
119       Crystal crystal, String filename) {
120     this.orix = orix;
121     this.oriy = oriy;
122     this.oriz = oriz;
123     this.extx = extx;
124     this.exty = exty;
125     this.extz = extz;
126     this.nx = nx;
127     this.ny = ny;
128     this.nz = nz;
129     this.crystal = crystal;
130     this.filename = filename;
131     this.stride = 2;
132   }
133 
134   /**
135    * Set the stepping across the array (e.g. 2 if data is separated by 1 space)
136    *
137    * @param stride the step size desired
138    */
139   public void setStride(int stride) {
140     this.stride = stride;
141   }
142 
143   /**
144    * Write data to file (does not normalize).
145    *
146    * @param data map data to write out
147    */
148   public void write(double[] data) {
149     write(data, false);
150   }
151 
152   /**
153    * write data to file, does not normalize
154    *
155    * @param data map data to write out
156    * @param norm should the data be normalized by mean/sd?
157    */
158   public void write(double[] data, boolean norm) {
159     ByteOrder b = ByteOrder.nativeOrder();
160     FileOutputStream fos;
161     DataOutputStream dos;
162 
163     double min = Double.POSITIVE_INFINITY;
164     double max = Double.NEGATIVE_INFINITY;
165     double mean = 0.0;
166     double sd = 0.0;
167 
168     int n = 0;
169     for (int k = 0; k < extz; k++) {
170       for (int j = 0; j < exty; j++) {
171         for (int i = 0; i < extx; i++) {
172           int index = stride * (i + extx * (j + exty * k));
173           n++;
174           if (data[index] < min) {
175             min = data[index];
176           }
177           if (data[index] > max) {
178             max = data[index];
179           }
180           mean += (data[index] - mean) / n;
181         }
182       }
183     }
184 
185     n = 0;
186     for (int k = 0; k < extz; k++) {
187       for (int j = 0; j < exty; j++) {
188         for (int i = 0; i < extx; i++) {
189           int index = stride * (i + extx * (j + exty * k));
190           sd += pow(data[index] - mean, 2.0);
191           n++;
192         }
193       }
194     }
195     sd = sqrt(sd / n);
196 
197     if (norm) {
198       for (int k = 0; k < extz; k++) {
199         for (int j = 0; j < exty; j++) {
200           for (int i = 0; i < extx; i++) {
201             int index = stride * (i + extx * (j + exty * k));
202             data[index] = (data[index] - mean) / sd;
203           }
204         }
205       }
206       // recurse
207       write(data, false);
208     }
209 
210     File file = version(new File(filename));
211     String name = file.getName();
212 
213     try {
214       if (logger.isLoggable(Level.INFO)) {
215         StringBuilder sb = new StringBuilder();
216         sb.append(format("\n Writing CCP4 map file: \"%s\"", name));
217         sb.append(format("\n  Map min: %g max: %g mean: %g standard dev.: %g", min, max, mean, sd));
218         logger.info(sb.toString());
219       }
220 
221       fos = new FileOutputStream(file);
222       dos = new DataOutputStream(fos);
223 
224       byte[] bytes = new byte[2048];
225       int offset = 0;
226 
227       int imapdata;
228       float fmapdata;
229       String mapstr;
230 
231       // header
232       ByteBuffer bb = ByteBuffer.wrap(bytes);
233       bb.order(b).putInt(extx);
234       bb.order(b).putInt(exty);
235       bb.order(b).putInt(extz);
236 
237       // mode (2 = reals, only one we accept)
238       bb.order(b).putInt(2);
239 
240       bb.order(b).putInt(orix);
241       bb.order(b).putInt(oriy);
242       bb.order(b).putInt(oriz);
243       bb.order(b).putInt(nx);
244       bb.order(b).putInt(ny);
245       bb.order(b).putInt(nz);
246 
247       bb.order(b).putFloat((float) crystal.a);
248       bb.order(b).putFloat((float) crystal.b);
249       bb.order(b).putFloat((float) crystal.c);
250       bb.order(b).putFloat((float) crystal.alpha);
251       bb.order(b).putFloat((float) crystal.beta);
252       bb.order(b).putFloat((float) crystal.gamma);
253 
254       bb.order(b).putInt(1);
255       bb.order(b).putInt(2);
256       bb.order(b).putInt(3);
257 
258       bb.order(b).putFloat((float) min);
259       bb.order(b).putFloat((float) max);
260       bb.order(b).putFloat((float) mean);
261 
262       bb.order(b).putInt(crystal.spaceGroup.number);
263       // bb.order(b).putInt(1);
264 
265       // symmetry bytes - should set this up at some point
266       // imapdata = swap ? ByteSwap.swap(320) : 320;
267       bb.order(b).putInt(80);
268 
269       bb.order(b).putInt(0);
270 
271       for (int i = 0; i < 12; i++) {
272         bb.order(b).putFloat(0.0f);
273       }
274 
275       for (int i = 0; i < 15; i++) {
276         bb.order(b).putInt(0);
277       }
278       dos.write(bytes, offset, 208);
279       bb.rewind();
280 
281       mapstr = "MAP ";
282       dos.writeBytes(mapstr);
283 
284       // machine code: double, float, int, uchar
285       // 0x4441 for LE, 0x1111 for BE
286       if (ByteOrder.nativeOrder().equals(ByteOrder.LITTLE_ENDIAN)) {
287         imapdata = 0x4441;
288       } else {
289         imapdata = 0x1111;
290       }
291       bb.order(b).putInt(imapdata);
292 
293       bb.order(b).putFloat((float) sd);
294 
295       bb.order(b).putInt(1);
296       dos.write(bytes, offset, 12);
297 
298       StringBuilder sb = new StringBuilder();
299       sb.append("map data from ffx");
300       while (sb.length() < 80) {
301         sb.append(" ");
302       }
303       dos.writeBytes(sb.toString());
304 
305       sb = new StringBuilder();
306       while (sb.length() < 80) {
307         sb.append(" ");
308       }
309       for (int i = 0; i < 9; i++) {
310         dos.writeBytes(sb.toString());
311       }
312 
313       sb = new StringBuilder();
314       sb.append("x,y,z");
315       while (sb.length() < 80) {
316         sb.append(" ");
317       }
318       dos.writeBytes(sb.toString());
319 
320       bb.rewind();
321       for (int k = 0; k < extz; k++) {
322         for (int j = 0; j < exty; j++) {
323           for (int i = 0; i < extx; i++) {
324             int index = stride * (i + extx * (j + exty * k));
325             fmapdata = (float) data[index];
326             bb.order(b).putFloat(fmapdata);
327             if (!bb.hasRemaining()) {
328               dos.write(bytes);
329               bb.rewind();
330             }
331           }
332         }
333       }
334       if (bb.position() > 0) {
335         dos.write(bytes);
336         bb.rewind();
337       }
338 
339       dos.close();
340     } catch (Exception e) {
341       String message = "Fatal exception evaluating structure factors.\n";
342       logger.log(Level.SEVERE, message, e);
343     }
344   }
345 }