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