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.SpaceGroupDefinitions;
45  import ffx.crystal.SpaceGroupInfo;
46  import ffx.xray.DiffractionRefinementData;
47  import org.apache.commons.configuration2.CompositeConfiguration;
48  
49  import java.io.BufferedReader;
50  import java.io.File;
51  import java.io.FileReader;
52  import java.io.IOException;
53  import java.util.logging.Level;
54  import java.util.logging.Logger;
55  
56  import static ffx.crystal.SpaceGroupInfo.spaceGroupNames;
57  import static java.lang.Double.parseDouble;
58  import static java.lang.Integer.parseInt;
59  import static java.lang.String.format;
60  import static org.apache.commons.math3.util.FastMath.min;
61  
62  /**
63   * CNSFilter class.
64   *
65   * @author Timothy D. Fenn
66   * @since 1.0
67   */
68  public class CNSFilter implements DiffractionFileFilter {
69  
70    private static final Logger logger = Logger.getLogger(CNSFilter.class.getName());
71  
72    private final double[] cell = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0};
73    private double resHigh = -1.0;
74    private String spaceGroupName = null;
75    private int spaceGroupNum = -1;
76  
77    /** Constructor for CNSFilter. */
78    public CNSFilter() {
79    }
80  
81    /** {@inheritDoc} */
82    @Override
83    public ReflectionList getReflectionList(File cnsFile) {
84      return getReflectionList(cnsFile, null);
85    }
86  
87    /** {@inheritDoc} */
88    @Override
89    public ReflectionList getReflectionList(File cnsFile, CompositeConfiguration properties) {
90      try (BufferedReader br = new BufferedReader(new FileReader(cnsFile))) {
91        String string;
92        while ((string = br.readLine()) != null) {
93          String[] strArray = string.split("\\s+");
94          if (strArray[0].equalsIgnoreCase("{")) {
95            if (strArray[1].toLowerCase().startsWith("sg=")) {
96              spaceGroupName = strArray[1].substring(3);
97              cell[0] = parseDouble(strArray[2].substring(2));
98              cell[1] = parseDouble(strArray[3].substring(2));
99              cell[2] = parseDouble(strArray[4].substring(2));
100             cell[3] = parseDouble(strArray[5].substring(6));
101             cell[4] = parseDouble(strArray[6].substring(5));
102             cell[5] = parseDouble(strArray[7].substring(6));
103           }
104         } else if (strArray[0].equalsIgnoreCase("CRYST1")) {
105           cell[0] = parseDouble(strArray[1]);
106           cell[1] = parseDouble(strArray[2]);
107           cell[2] = parseDouble(strArray[3]);
108           cell[3] = parseDouble(strArray[4]);
109           cell[4] = parseDouble(strArray[5]);
110           cell[5] = parseDouble(strArray[6]);
111           spaceGroupName = SpaceGroupInfo.pdb2ShortName(string.substring(55, 65));
112         } else if (strArray[0].toLowerCase().startsWith("inde")) {
113           break;
114         }
115       }
116     } catch (IOException e) {
117       String message = " CNS IO Exception.";
118       logger.log(Level.WARNING, message, e);
119       return null;
120     }
121 
122     Resolution resolution = null;
123     if (properties != null) {
124       resolution = Resolution.checkProperties(properties);
125       resHigh = resolution.resolution;
126     }
127 
128     if (spaceGroupName != null) {
129       spaceGroupNum = SpaceGroupDefinitions.spaceGroupNumber(spaceGroupName);
130     }
131 
132     if (spaceGroupNum < 0 || cell[0] < 0 || resolution == null) {
133       logger.info(
134           " The CNS file contains insufficient information to generate the reflection list.");
135       return null;
136     }
137 
138     if (logger.isLoggable(Level.INFO)) {
139       StringBuilder sb = new StringBuilder();
140       sb.append(format("\n Opening %s\n", cnsFile.getName()));
141       sb.append(" Setting up Reflection List based on CNS:\n");
142       sb.append(format("  Spacegroup #: %d (name: %s)\n",
143           spaceGroupNum, spaceGroupNames[spaceGroupNum - 1]));
144       sb.append(format("  Resolution:   %8.3f\n", resHigh));
145       sb.append(format("  Cell:         %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
146           cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]));
147       logger.info(sb.toString());
148     }
149 
150     Crystal crystal = new Crystal(cell[0], cell[1], cell[2],
151         cell[3], cell[4], cell[5], spaceGroupNames[spaceGroupNum - 1]);
152 
153     return new ReflectionList(crystal, resolution, properties);
154   }
155 
156   /** {@inheritDoc} */
157   @Override
158   public double getResolution(File cnsFile, Crystal crystal) {
159     double res = Double.POSITIVE_INFINITY;
160 
161     try (BufferedReader br = new BufferedReader(new FileReader(cnsFile))) {
162       HKL hkl = new HKL();
163       String string;
164       while ((string = br.readLine()) != null) {
165         String[] strArray = string.split("\\s+");
166         for (int i = 0; i < strArray.length; i++) {
167           if (strArray[i].toLowerCase().startsWith("inde")) {
168             if (i < strArray.length - 3) {
169               int ih = parseInt(strArray[i + 1]);
170               int ik = parseInt(strArray[i + 2]);
171               int il = parseInt(strArray[i + 3]);
172               hkl.setH(ih);
173               hkl.setK(ik);
174               hkl.setL(il);
175               res = min(res, crystal.res(hkl));
176             }
177           }
178         }
179       }
180     } catch (IOException e) {
181       String message = " CNS IO Exception.";
182       logger.log(Level.WARNING, message, e);
183       return -1.0;
184     }
185 
186     return res;
187   }
188 
189   /** {@inheritDoc} */
190   @Override
191   public boolean readFile(
192       File cnsFile,
193       ReflectionList reflectionList,
194       DiffractionRefinementData refinementData,
195       CompositeConfiguration properties) {
196 
197     int nRead, nRes, nIgnore, nFriedel, nCut;
198     boolean transpose = false;
199 
200     StringBuilder sb = new StringBuilder();
201     sb.append(format("\n Opening %s\n", cnsFile.getName()));
202 
203     if (refinementData.rFreeFlag < 0) {
204       refinementData.setFreeRFlag(1);
205       sb.append(format(" Setting R free flag to CNS default: %d\n", refinementData.rFreeFlag));
206     }
207 
208     // column identifiers
209     String foString = null;
210     String sigFoString = null;
211     String rFreeString = null;
212     int ih, ik, il, free;
213     double fo, sigFo;
214 
215     try (BufferedReader br = new BufferedReader(new FileReader(cnsFile))) {
216       ih = ik = il = free = -1;
217       fo = sigFo = -1.0;
218 
219       // Check if HKLs need to be transposed or not.
220       HKL mate = new HKL();
221       int nPosIgnore = 0;
222       int nTransIgnore = 0;
223 
224       String string;
225       while ((string = br.readLine()) != null) {
226         String[] strArray = string.split("\\s+");
227 
228         for (int i = 0; i < strArray.length; i++) {
229           if (strArray[i].toLowerCase().startsWith("inde")) {
230             if (i < strArray.length - 3) {
231               ih = parseInt(strArray[i + 1]);
232               ik = parseInt(strArray[i + 2]);
233               il = parseInt(strArray[i + 3]);
234               reflectionList.findSymHKL(ih, ik, il, mate, false);
235               HKL hklpos = reflectionList.getHKL(mate);
236               if (hklpos == null) {
237                 nPosIgnore++;
238               }
239 
240               reflectionList.findSymHKL(ih, ik, il, mate, true);
241               HKL hkltrans = reflectionList.getHKL(mate);
242               if (hkltrans == null) {
243                 nTransIgnore++;
244               }
245             }
246           }
247         }
248       }
249       if (nPosIgnore > nTransIgnore) {
250         transpose = true;
251       }
252 
253       if (properties != null) {
254         foString = properties.getString("fostring", null);
255         sigFoString = properties.getString("sigfostring", null);
256         rFreeString = properties.getString("rfreestring", null);
257       }
258     } catch (IOException e) {
259       String message = "CNS IO Exception.";
260       logger.log(Level.WARNING, message, e);
261       return false;
262     }
263 
264     boolean hasHKL, hasFo, hasSigFo, hasFree;
265     hasHKL = hasFo = hasSigFo = hasFree = false;
266 
267     // reopen to start at beginning
268     try (BufferedReader br = new BufferedReader(new FileReader(cnsFile))) {
269 
270       // read in data
271       double[][] anofSigF = new double[refinementData.n][4];
272       for (int i = 0; i < refinementData.n; i++) {
273         anofSigF[i][0] = anofSigF[i][1] = anofSigF[i][2] = anofSigF[i][3] = Double.NaN;
274       }
275       nRead = nRes = nIgnore = nFriedel = nCut = 0;
276 
277       String string;
278       HKL mate = new HKL();
279       while ((string = br.readLine()) != null) {
280         String[] strArray = string.split("\\s+");
281         for (int i = 0; i < strArray.length; i++) {
282           if (strArray[i].toLowerCase().startsWith("inde")) {
283             if (hasHKL && hasFo && hasSigFo && hasFree) {
284               boolean friedel = reflectionList.findSymHKL(ih, ik, il, mate, transpose);
285               HKL hkl = reflectionList.getHKL(mate);
286               if (hkl != null) {
287                 if (refinementData.fSigFCutoff > 0.0 && (fo / sigFo) < refinementData.fSigFCutoff) {
288                   nCut++;
289                 } else if (friedel) {
290                   anofSigF[hkl.getIndex()][2] = fo;
291                   anofSigF[hkl.getIndex()][3] = sigFo;
292                   nFriedel++;
293                 } else {
294                   anofSigF[hkl.getIndex()][0] = fo;
295                   anofSigF[hkl.getIndex()][1] = sigFo;
296                 }
297                 refinementData.setFreeR(hkl.getIndex(), free);
298                 nRead++;
299               } else {
300                 HKL tmp = new HKL(ih, ik, il);
301                 if (!reflectionList.resolution.inInverseResSqRange(
302                     reflectionList.crystal.invressq(tmp))) {
303                   nRes++;
304                 } else {
305                   nIgnore++;
306                 }
307               }
308             }
309             hasHKL = false;
310             hasFo = false;
311             hasSigFo = false;
312             hasFree = false;
313             if (i < strArray.length - 3) {
314               ih = parseInt(strArray[i + 1]);
315               ik = parseInt(strArray[i + 2]);
316               il = parseInt(strArray[i + 3]);
317               hasHKL = true;
318             }
319           }
320           if (strArray[i].toLowerCase().startsWith("fobs=")
321               || strArray[i].equalsIgnoreCase(foString + "=")) {
322             fo = parseDouble(strArray[i + 1]);
323             hasFo = true;
324           }
325           if (strArray[i].toLowerCase().startsWith("sigma=")
326               || strArray[i].equalsIgnoreCase(sigFoString + "=")) {
327             sigFo = parseDouble(strArray[i + 1]);
328             hasSigFo = true;
329           }
330           if (strArray[i].toLowerCase().startsWith("test=")
331               || strArray[i].equalsIgnoreCase(rFreeString + "=")) {
332             free = parseInt(strArray[i + 1]);
333             hasFree = true;
334           }
335         }
336       }
337       // Set up fsigf from F+ and F-.
338       refinementData.generateFsigFfromAnomalousFsigF(anofSigF);
339     } catch (IOException e) {
340       String message = "CNS IO Exception.";
341       logger.log(Level.WARNING, message, e);
342       return false;
343     }
344 
345     if (logger.isLoggable(Level.INFO)) {
346       sb.append(format(" HKL read in:                             %d\n", nRead));
347       sb.append(format(" HKL read as friedel mates:               %d\n", nFriedel));
348       sb.append(format(" HKL not read in (too high resolution):   %d\n", nRes));
349       sb.append(format(" HKL not read in (not in internal list?): %d\n", nIgnore));
350       sb.append(format(" HKL not read in (F/sigF cutoff):         %d\n", nCut));
351       sb.append(
352           format(" HKL in internal list:                    %d\n", reflectionList.hklList.size()));
353       logger.info(sb.toString());
354     }
355 
356     return true;
357   }
358 }