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.crystal;
39  
40  import org.apache.commons.configuration2.CompositeConfiguration;
41  
42  import javax.annotation.Nullable;
43  import java.util.ArrayList;
44  import java.util.HashMap;
45  import java.util.Map.Entry;
46  import java.util.logging.Logger;
47  
48  import static ffx.crystal.SymOp.applySymRot;
49  import static ffx.crystal.SymOp.applyTransSymRot;
50  import static ffx.numerics.math.ScalarMath.mod;
51  import static org.apache.commons.math3.util.FastMath.PI;
52  import static org.apache.commons.math3.util.FastMath.cos;
53  import static org.apache.commons.math3.util.FastMath.floor;
54  import static org.apache.commons.math3.util.FastMath.max;
55  import static org.apache.commons.math3.util.FastMath.min;
56  import static org.apache.commons.math3.util.FastMath.rint;
57  
58  /**
59   * Class to represent a reflection list.
60   *
61   * @author Timothy D. Fenn
62   * @see <a href="http://dx.doi.org/10.1107/S0021889802013420" target="_blank"> Cowtan, K. 2002.
63   * Generic representation and evaluation of properties as a function of position in reciprocal
64   * space. J. Appl. Cryst. 35:655-663. </a>
65   * @since 1.0
66   */
67  public class ReflectionList {
68  
69    private static final Logger logger = Logger.getLogger(ReflectionList.class.getName());
70  
71    /**
72     * The HKL list.
73     */
74    public final ArrayList<HKL> hklList = new ArrayList<>();
75    /**
76     * The Crystal instance.
77     */
78    public final Crystal crystal;
79    /**
80     * The space group.
81     */
82    public final SpaceGroup spaceGroup;
83    /**
84     * Resolution instance.
85     */
86    public final Resolution resolution;
87    /**
88     * String to HKL look-up.
89     */
90    final HashMap<String, HKL> hklMap = new HashMap<>();
91    /**
92     * For binning reflections based on resolution
93     */
94    public int nBins = 10;
95    /**
96     * Histogram.
97     */
98    private final double[] histogram = new double[1001];
99    /**
100    * Minimum resolution.
101    */
102   private double minResolution;
103   /**
104    * Maximum resolution.
105    */
106   private double maxResolution;
107 
108   /**
109    * Constructor for ReflectionList.
110    *
111    * @param crystal    a {@link ffx.crystal.Crystal} object.
112    * @param resolution a {@link ffx.crystal.Resolution} object.
113    */
114   public ReflectionList(Crystal crystal, Resolution resolution) {
115     this(crystal, resolution, null);
116   }
117 
118   /**
119    * Constructor for ReflectionList.
120    *
121    * @param crystal    a {@link ffx.crystal.Crystal} object.
122    * @param resolution a {@link ffx.crystal.Resolution} object.
123    * @param properties a {@link org.apache.commons.configuration2.CompositeConfiguration} object.
124    */
125   public ReflectionList(Crystal crystal, Resolution resolution, CompositeConfiguration properties) {
126     this.crystal = crystal;
127     this.spaceGroup = crystal.spaceGroup;
128     this.resolution = resolution;
129 
130     int hMax = (int) (this.crystal.a / this.resolution.resolutionLimit());
131     int kMax = (int) (this.crystal.b / this.resolution.resolutionLimit());
132     int lMax = (int) (this.crystal.c / this.resolution.resolutionLimit());
133 
134     minResolution = Double.POSITIVE_INFINITY;
135     maxResolution = Double.NEGATIVE_INFINITY;
136     int n = 0;
137 
138     HKL hkl = new HKL();
139     for (int h = -hMax; h <= hMax; h++) {
140       hkl.setH(h);
141       for (int k = -kMax; k <= kMax; k++) {
142         hkl.setK(k);
143         for (int l = -lMax; l <= lMax; l++) {
144           hkl.setL(l);
145 
146           double res = this.crystal.invressq(hkl);
147           getEpsilon(hkl);
148           LaueSystem laueSystem = spaceGroup.laueSystem;
149           if (laueSystem.checkRestrictions(h, k, l)
150               && resolution.inInverseResSqRange(res)
151               && !hkl.sysAbs()) {
152             minResolution = min(res, minResolution);
153             maxResolution = max(res, maxResolution);
154             String s = (h + "_" + k + "_" + l).intern();
155             hklMap.put(s, new HKL(hkl.getH(), hkl.getK(), hkl.getL(), hkl.getEpsilon(), hkl.allowed));
156             n++;
157           }
158         }
159       }
160     }
161 
162     n = 0;
163     for (Entry<String, HKL> entry : hklMap.entrySet()) {
164       HKL ih = entry.getValue();
165       ih.setIndex(n);
166       hklList.add(ih);
167       n++;
168     }
169 
170     // Set up the resolution bins first build a histogram.
171     for (HKL ih : hklList) {
172       double r = (this.crystal.invressq(ih) - minResolution) / (maxResolution - minResolution);
173       int i = (int) (min(r, 0.999) * 1000.0);
174       histogram[i + 1] += 1.0;
175     }
176 
177     // Convert to cumulative histogram
178     for (int i = 1; i < histogram.length; i++) {
179       histogram[i] += histogram[i - 1];
180     }
181     for (int i = 0; i < histogram.length; i++) {
182       histogram[i] /= histogram[histogram.length - 1];
183     }
184 
185     // Assign each reflection to a bin in the range (0-nbins)
186     setResolutionBins(properties);
187   }
188 
189   /**
190    * Constructor for ReflectionList.
191    *
192    * @param a          a double.
193    * @param b          a double.
194    * @param c          a double.
195    * @param alpha      a double.
196    * @param beta       a double.
197    * @param gamma      a double.
198    * @param sg         a {@link java.lang.String} object.
199    * @param resolution a double.
200    */
201   public ReflectionList(double a, double b, double c, double alpha, double beta, double gamma,
202                         String sg, double resolution) {
203     this(new Crystal(a, b, c, alpha, beta, gamma, sg), new Resolution(resolution));
204   }
205 
206   /**
207    * findSymHKL
208    *
209    * @param hkl  a {@link ffx.crystal.HKL} object.
210    * @param mate a {@link ffx.crystal.HKL} object.
211    * @return a boolean.
212    */
213   public boolean findSymHKL(HKL hkl, HKL mate) {
214     return findSymHKL(hkl, mate, false);
215   }
216 
217   /**
218    * findSymHKL
219    *
220    * @param h    an int.
221    * @param k    an int.
222    * @param l    an int.
223    * @param mate a {@link ffx.crystal.HKL} object.
224    * @return a boolean.
225    */
226   public boolean findSymHKL(int h, int k, int l, HKL mate) {
227     return findSymHKL(new HKL(h, k, l), mate, false);
228   }
229 
230   /**
231    * findSymHKL
232    *
233    * @param h         an int.
234    * @param k         an int.
235    * @param l         an int.
236    * @param mate      a {@link ffx.crystal.HKL} object.
237    * @param transpose a boolean.
238    * @return a boolean.
239    */
240   public boolean findSymHKL(int h, int k, int l, HKL mate, boolean transpose) {
241     return findSymHKL(new HKL(h, k, l), mate, transpose);
242   }
243 
244   /**
245    * getHKL
246    *
247    * @param h an int.
248    * @param k an int.
249    * @param l an int.
250    * @return a {@link ffx.crystal.HKL} object.
251    */
252   public HKL getHKL(int h, int k, int l) {
253     String s = (h + "_" + k + "_" + l);
254     return hklMap.get(s);
255   }
256 
257   /**
258    * getHKL
259    *
260    * @param hkl a {@link ffx.crystal.HKL} object.
261    * @return a {@link ffx.crystal.HKL} object.
262    */
263   public HKL getHKL(HKL hkl) {
264     return getHKL(hkl.getH(), hkl.getK(), hkl.getL());
265   }
266 
267   /**
268    * Get the maximum resolution.
269    *
270    * @return Maximum resolution
271    */
272   public double getMaxResolution() {
273     return maxResolution;
274   }
275 
276   /**
277    * Get the minimum resolution.
278    *
279    * @return Minimum resolution
280    */
281   public double getMinResolution() {
282     return minResolution;
283   }
284 
285   /**
286    * hasHKL
287    *
288    * @param hkl a {@link ffx.crystal.HKL} object.
289    * @return a boolean.
290    */
291   public boolean hasHKL(HKL hkl) {
292     return hasHKL(hkl.getH(), hkl.getK(), hkl.getL());
293   }
294 
295   /**
296    * ordinal
297    *
298    * @param s a double.
299    * @return a double.
300    */
301   public final double ordinal(double s) {
302     double r = (s - minResolution) / (maxResolution - minResolution);
303     r = min(r, 0.999) * 1000.0;
304     int i = (int) r;
305     r -= floor(r);
306     return ((1.0 - r) * histogram[i] + r * histogram[i + 1]);
307   }
308 
309   /**
310    * {@inheritDoc}
311    */
312   @Override
313   public String toString() {
314     return " Reflection list with "
315         + this.hklList.size()
316         + " reflections, spacegroup "
317         + this.spaceGroup.shortName
318         + " resolution limit: "
319         + resolution.resolutionLimit();
320   }
321 
322   /**
323    * findSymHKL
324    *
325    * @param hkl       a {@link ffx.crystal.HKL} object.
326    * @param mate      a {@link ffx.crystal.HKL} object.
327    * @param transpose a boolean.
328    * @return a boolean.
329    */
330   private boolean findSymHKL(HKL hkl, HKL mate, boolean transpose) {
331     int nsym = spaceGroup.numPrimitiveSymEquiv;
332 
333     for (int i = 0; i < nsym; i++) {
334       if (transpose) {
335         applyTransSymRot(hkl, mate, spaceGroup.symOps.get(i));
336       } else {
337         applySymRot(hkl, mate, spaceGroup.symOps.get(i));
338       }
339 
340       LaueSystem laueSystem = spaceGroup.laueSystem;
341       if (laueSystem.checkRestrictions(mate.getH(), mate.getK(), mate.getL())) {
342         return false;
343       }
344       if (laueSystem.checkRestrictions(-mate.getH(), -mate.getK(), -mate.getL())) {
345         mate.setH(-mate.getH());
346         mate.setK(-mate.getK());
347         mate.setL(-mate.getL());
348         return true;
349       }
350     }
351 
352     mate.setH(hkl.getH());
353     mate.setK(hkl.getK());
354     mate.setL(hkl.getL());
355     return false;
356   }
357 
358   /**
359    * hasHKL
360    *
361    * @param h an int.
362    * @param k an int.
363    * @param l an int.
364    * @return a boolean.
365    */
366   private boolean hasHKL(int h, int k, int l) {
367     String s = (h + "_" + k + "_" + l);
368     return hklMap.containsKey(s);
369   }
370 
371   private void getEpsilon(HKL hkl) {
372     int epsilon = 1;
373     int allowed = 255;
374 
375     int nSym = spaceGroup.symOps.size();
376     for (int i = 1; i < nSym; i++) {
377       HKL mate = new HKL();
378       SymOp symOp = spaceGroup.symOps.get(i);
379       applySymRot(hkl, mate, symOp);
380       double shift = symOp.symPhaseShift(hkl);
381       if (mate.equals(hkl)) {
382         if (cos(shift) > 0.999) {
383           epsilon++;
384         } else {
385           allowed = 0;
386           epsilon = 0;
387           break;
388         }
389       } else if (mate.equals(hkl.neg())) {
390         // centric reflection
391         allowed = (int) rint(mod(-0.5 * shift, PI) / (PI / HKL.ndiv));
392       }
393     }
394     if (hkl.getH() == 0 && hkl.getK() == 0 && hkl.getL() == 0) {
395       allowed = 0;
396     }
397 
398     hkl.setEpsilon(epsilon);
399     hkl.setAllowed(allowed);
400   }
401 
402   private void setResolutionBins(@Nullable CompositeConfiguration properties) {
403     if (properties != null) {
404       nBins = properties.getInt("reflection-bins", 10);
405     }
406     double nBinsDouble = nBins;
407     for (HKL ih : hklList) {
408       int bin = (int) (nBinsDouble * ordinal(crystal.invressq(ih)));
409       ih.bin = min(bin, nBins - 1);
410     }
411   }
412 }