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