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;
39  
40  import static java.lang.Double.isNaN;
41  import static java.lang.String.format;
42  import static org.apache.commons.math3.util.FastMath.PI;
43  import static org.apache.commons.math3.util.FastMath.sqrt;
44  
45  import ffx.crystal.ReflectionList;
46  import ffx.numerics.math.ComplexNumber;
47  import java.util.Random;
48  import java.util.logging.Level;
49  import java.util.logging.Logger;
50  import org.apache.commons.configuration2.CompositeConfiguration;
51  
52  /**
53   * DiffractionRefinementData class.
54   *
55   * @author Timothy D. Fenn
56   * @since 1.0
57   */
58  public class DiffractionRefinementData {
59  
60    private static final Logger logger = Logger.getLogger(DiffractionRefinementData.class.getName());
61    /** Number of reflections in the data set. */
62    public final int n;
63    /** Array of R free flags; */
64    public final int[] freeR;
65    /** Calculated atomic structure factors. */
66    public final double[][] fc;
67    /** Calculated bulk solvent structure factors. */
68    public final double[][] fs;
69    /** Figure of merit and phase. */
70    public final double[][] fomPhi;
71    /** mFo - DFc coefficients. */
72    public final double[][] foFc1;
73  
74    public final double fSigFCutoff;
75    /** Number of scale parameters. */
76    final int nScale;
77    /** 2D array of F/sigF data. */
78    final double[][] fSigF;
79    /** Scaled sum of Fc and Fs */
80    final double[][] fcTot;
81    /** 2mFo - DFc coefficients. */
82    final double[][] foFc2;
83    /** Derivatives with respect to Fc. */
84    final double[][] dFc;
85    /** Derivatives with respect to Fs. */
86    final double[][] dFs;
87    /** Number of resolution bins. */
88    final int nBins;
89    /** Spine scaling coefficients. */
90    public double[] spline;
91    /** SigmaA coefficient - s. */
92    public double[] sigmaA;
93    /** SigmaA coefficient - w. */
94    public double[] sigmaW;
95    /**
96     * Duplicated settings - these are also in DiffractionData, but duplicated here until settings are
97     * put in their own class.
98     */
99    public int rFreeFlag;
100   /** Log Likelihoods. */
101   double llkR, llkF;
102   /** Reciprocal space reference for structure factor calculations and computing derivatives. */
103   CrystalReciprocalSpace crystalReciprocalSpaceFc;
104   /**
105    * Reciprocal space reference for bulk solvent structure factor calculations and computing
106    * derivatives.
107    */
108   CrystalReciprocalSpace crystalReciprocalSpaceFs;
109   /** Scaled, E-like Fc. */
110   double[] esqFc;
111   /** Scaled, E-like Fo. */
112   double[] esqFo;
113   /** Bulk solvent parameters. */
114   double solventA, solventB;
115   /** bulk solvent scale and Ueq */
116   double bulkSolventK, bulkSolventUeq;
117   /** Overall model scale. */
118   double modelScaleK;
119   /** model anisotropic B */
120   double[] modelAnisoB = new double[6];
121 
122   /**
123    * allocate data given a {@link ffx.crystal.ReflectionList}
124    *
125    * @param properties configuration properties
126    * @param reflectionList {@link ffx.crystal.ReflectionList} to use to allocate data
127    */
128   public DiffractionRefinementData(
129       CompositeConfiguration properties, ReflectionList reflectionList) {
130 
131     int rflag = -1;
132     if (properties != null) {
133       rflag = properties.getInt("rfree-flag", -1);
134       fSigFCutoff = properties.getDouble("f-sigf-cutoff", -1.0);
135     } else {
136       fSigFCutoff = -1.0;
137     }
138 
139     this.n = reflectionList.hklList.size();
140     this.nScale = reflectionList.crystal.scaleN;
141     this.rFreeFlag = rflag;
142     this.nBins = reflectionList.nBins;
143     fSigF = new double[n][2];
144     freeR = new int[n];
145     fc = new double[n][2];
146     fs = new double[n][2];
147     fcTot = new double[n][2];
148     fomPhi = new double[n][2];
149     foFc2 = new double[n][2];
150     foFc1 = new double[n][2];
151     dFc = new double[n][2];
152     dFs = new double[n][2];
153 
154     for (int i = 0; i < n; i++) {
155       fSigF[i][0] = fSigF[i][1] = Double.NaN;
156       fcTot[i][0] = fcTot[i][1] = Double.NaN;
157     }
158 
159     spline = new double[nBins * 2];
160     sigmaA = new double[nBins];
161     sigmaW = new double[nBins];
162     esqFc = new double[nBins];
163     esqFo = new double[nBins];
164     for (int i = 0; i < nBins; i++) {
165       spline[i] = spline[i + nBins] = sigmaA[i] = esqFc[i] = esqFo[i] = 1.0;
166       sigmaW[i] = 0.05;
167     }
168 
169     // Initial guess for scaling/bulk solvent correction.
170     bulkSolventK = 0.33;
171     bulkSolventUeq = 50.0 / (8.0 * PI * PI);
172     modelScaleK = 0.0;
173   }
174 
175   /**
176    * FoFc2F
177    *
178    * @param i a int.
179    * @return a double.
180    */
181   public double FoFc2F(int i) {
182     ComplexNumber c = new ComplexNumber(foFc2[i][0], foFc2[i][1]);
183     return c.abs();
184   }
185 
186   /**
187    * FoFc2Phi
188    *
189    * @param i a int.
190    * @return a double.
191    */
192   public double FoFc2Phi(int i) {
193     ComplexNumber c = new ComplexNumber(foFc2[i][0], foFc2[i][1]);
194     return c.phase();
195   }
196 
197   /**
198    * get the amplitude of a complex Fc
199    *
200    * @param i reflection to get
201    * @return amplitude of Fc
202    */
203   public double fcF(int i) {
204     ComplexNumber c = new ComplexNumber(fc[i][0], fc[i][1]);
205     return c.abs();
206   }
207 
208   /**
209    * get the phase of a complex Fc
210    *
211    * @param i reflection to get
212    * @return phase of Fc
213    */
214   public double fcPhi(int i) {
215     ComplexNumber c = new ComplexNumber(fc[i][0], fc[i][1]);
216     return c.phase();
217   }
218 
219   /**
220    * fcTotF
221    *
222    * @param i a int.
223    * @return a double.
224    */
225   public double fcTotF(int i) {
226     ComplexNumber c = new ComplexNumber(fcTot[i][0], fcTot[i][1]);
227     return c.abs();
228   }
229 
230   /**
231    * fcTotPhi
232    *
233    * @param i a int.
234    * @return a double.
235    */
236   public double fcTotPhi(int i) {
237     ComplexNumber c = new ComplexNumber(fcTot[i][0], fcTot[i][1]);
238     return c.phase();
239   }
240 
241   /**
242    * foFc1F
243    *
244    * @param i a int.
245    * @return a double.
246    */
247   public double foFc1F(int i) {
248     ComplexNumber c = new ComplexNumber(foFc1[i][0], foFc1[i][1]);
249     return c.abs();
250   }
251 
252   /**
253    * foFc1Phi
254    *
255    * @param i a int.
256    * @return a double.
257    */
258   public double foFc1Phi(int i) {
259     ComplexNumber c = new ComplexNumber(foFc1[i][0], foFc1[i][1]);
260     return c.phase();
261   }
262 
263   /**
264    * fsF
265    *
266    * @param i a int.
267    * @return a double.
268    */
269   public double fsF(int i) {
270     ComplexNumber c = new ComplexNumber(fs[i][0], fs[i][1]);
271     return c.abs();
272   }
273 
274   /**
275    * fsPhi
276    *
277    * @param i a int.
278    * @return a double.
279    */
280   public double fsPhi(int i) {
281     ComplexNumber c = new ComplexNumber(fs[i][0], fs[i][1]);
282     return c.phase();
283   }
284 
285   /**
286    * Generate average F/sigF from anomalous F/sigF.
287    *
288    * @param anomalousFsigF an array of double.
289    */
290   public void generateFsigFfromAnomalousFsigF(double[][] anomalousFsigF) {
291     for (int i = 0; i < n; i++) {
292       if (isNaN(anomalousFsigF[i][0]) && isNaN(anomalousFsigF[i][2])) {
293         // Do Nothing.
294       } else if (isNaN(anomalousFsigF[i][0])) {
295         fSigF[i][0] = anomalousFsigF[i][2];
296         fSigF[i][1] = anomalousFsigF[i][3];
297       } else if (isNaN(anomalousFsigF[i][2])) {
298         fSigF[i][0] = anomalousFsigF[i][0];
299         fSigF[i][1] = anomalousFsigF[i][1];
300       } else {
301         fSigF[i][0] = (anomalousFsigF[i][0] + anomalousFsigF[i][2]) / 2.0;
302         fSigF[i][1] =
303             sqrt(
304                 anomalousFsigF[i][1] * anomalousFsigF[i][1]
305                     + anomalousFsigF[i][3] * anomalousFsigF[i][3]);
306       }
307     }
308   }
309 
310   /** Mark 5% of reflections for cross validation (R free flags). */
311   public void generateRFree() {
312     Random generator = new Random();
313     int free;
314     int nonfree;
315     int nfree = 0;
316 
317     if (rFreeFlag == 0) {
318       free = 0;
319       nonfree = 1;
320     } else {
321       free = 1;
322       nonfree = 0;
323     }
324 
325     for (int i = 0; i < n; i++) {
326       if (isNaN(fSigF[i][0])) {
327         freeR[i] = nonfree;
328         continue;
329       }
330 
331       int randomi = generator.nextInt(100);
332       if (randomi < 5) {
333         freeR[i] = free;
334         nfree++;
335       } else {
336         freeR[i] = nonfree;
337       }
338     }
339 
340     if (logger.isLoggable(Level.INFO)) {
341       logger.info(format(" Assigned %d of %d reflections to the R free set (5%%).\n", nfree, n));
342     }
343   }
344 
345   /**
346    * getF
347    *
348    * @param i a int.
349    * @return a double.
350    */
351   public double getF(int i) {
352     return fSigF[i][0];
353   }
354 
355   /**
356    * getFSigF
357    *
358    * @param i a int.
359    * @return an array of double.
360    */
361   public double[] getFSigF(int i) {
362     return fSigF[i];
363   }
364 
365   /**
366    * getFcTot
367    *
368    * @param i a int.
369    * @return a {@link ComplexNumber} object.
370    */
371   public ComplexNumber getFcTot(int i) {
372     return new ComplexNumber(fcTot[i][0], fcTot[i][1]);
373   }
374 
375   /**
376    * getFoFc1
377    *
378    * @param i a int.
379    * @return a {@link ComplexNumber} object.
380    */
381   public ComplexNumber getFoFc1(int i) {
382     return new ComplexNumber(foFc1[i][0], foFc1[i][1]);
383   }
384 
385   /**
386    * getFoFc1IP
387    *
388    * @param i a int.
389    * @param c a {@link ComplexNumber} object.
390    */
391   public void getFoFc1IP(int i, ComplexNumber c) {
392     c.re(foFc1[i][0]);
393     c.im(foFc1[i][1]);
394   }
395 
396   /**
397    * getFoFc2
398    *
399    * @param i a int.
400    * @return a {@link ComplexNumber} object.
401    */
402   public ComplexNumber getFoFc2(int i) {
403     return new ComplexNumber(foFc2[i][0], foFc2[i][1]);
404   }
405 
406   /**
407    * getFoFc2IP
408    *
409    * @param i a int.
410    * @param c a {@link ComplexNumber} object.
411    */
412   public void getFoFc2IP(int i, ComplexNumber c) {
413     c.re(foFc2[i][0]);
414     c.im(foFc2[i][1]);
415   }
416 
417   /**
418    * getFreeR
419    *
420    * @param i a int.
421    * @return a int.
422    */
423   public int getFreeR(int i) {
424     return freeR[i];
425   }
426 
427   /**
428    * getFs
429    *
430    * @param i a int.
431    * @return a {@link ComplexNumber} object.
432    */
433   public ComplexNumber getFs(int i) {
434     return new ComplexNumber(fs[i][0], fs[i][1]);
435   }
436 
437   /**
438    * getSigF
439    *
440    * @param i a int.
441    * @return a double.
442    */
443   public double getSigF(int i) {
444     return fSigF[i][1];
445   }
446 
447   /**
448    * Generate amplitudes from intensities.
449    *
450    * <p>This does NOT use French and Wilson caling, but just a simple square root.
451    */
452   public void intensitiesToAmplitudes() {
453     double tmp;
454 
455     for (int i = 0; i < n; i++) {
456       if (fSigF[i][0] > 0.0) {
457         tmp = fSigF[i][0];
458         fSigF[i][0] = sqrt(tmp);
459         if (fSigF[i][1] < tmp) {
460           fSigF[i][1] = fSigF[i][0] - sqrt(tmp - fSigF[i][1]);
461         } else {
462           fSigF[i][1] = fSigF[i][0];
463         }
464       } else if (!isNaN(fSigF[i][0])) {
465         fSigF[i][0] = 0.0;
466         fSigF[i][1] = 0.0;
467       }
468     }
469 
470     if (logger.isLoggable(Level.WARNING)) {
471       StringBuilder sb = new StringBuilder();
472       sb.append("\n Internally converting intensities to amplitudes.\n");
473       sb.append(" This does not use French & Wilson scaling.\n");
474       logger.info(sb.toString());
475     }
476   }
477 
478   /**
479    * isFreeR
480    *
481    * @param i a int.
482    * @param f a int.
483    * @return a boolean.
484    */
485   public boolean isFreeR(int i, int f) {
486     return (freeR[i] == f);
487   }
488 
489   /**
490    * return the current likelihood
491    *
492    * @return the free likelihood (Rfree based)
493    */
494   public double likelihoodFree() {
495     return llkF;
496   }
497 
498   /**
499    * return the current likelihood
500    *
501    * @return the work likelihood (non-Rfree based)
502    */
503   public double likelihoodWork() {
504     return llkR;
505   }
506 
507   /**
508    * Set amplitude (F).
509    *
510    * @param i reflection to set
511    * @param f value of F desired
512    */
513   public void setF(int i, double f) {
514     fSigF[i][0] = f;
515   }
516 
517   /**
518    * Set amplitude and sigF.
519    *
520    * @param i reflection to set
521    * @param f value of F desired
522    * @param sigf value of sigF desired
523    */
524   public void setFSigF(int i, double f, double sigf) {
525     fSigF[i][0] = f;
526     fSigF[i][1] = sigf;
527   }
528 
529   /**
530    * set complex Fc
531    *
532    * @param i reflection to set
533    * @param c {@link ComplexNumber} to set reflection to
534    */
535   public void setFc(int i, ComplexNumber c) {
536     fc[i][0] = c.re();
537     fc[i][1] = c.im();
538   }
539 
540   /**
541    * setFcTot
542    *
543    * @param i a int.
544    * @param c a {@link ComplexNumber} object.
545    */
546   public void setFcTot(int i, ComplexNumber c) {
547     fcTot[i][0] = c.re();
548     fcTot[i][1] = c.im();
549   }
550 
551   /**
552    * setFoFc1
553    *
554    * @param i a int.
555    * @param c a {@link ComplexNumber} object.
556    */
557   public void setFoFc1(int i, ComplexNumber c) {
558     foFc1[i][0] = c.re();
559     foFc1[i][1] = c.im();
560   }
561 
562   /**
563    * setFoFc2
564    *
565    * @param i a int.
566    * @param c a {@link ComplexNumber} object.
567    */
568   public void setFoFc2(int i, ComplexNumber c) {
569     foFc2[i][0] = c.re();
570     foFc2[i][1] = c.im();
571   }
572 
573   /**
574    * Set FreeR value flag of a reflection.
575    *
576    * @param i reflection to set
577    * @param f FreeR value to set reflection to
578    */
579   public void setFreeR(int i, int f) {
580     freeR[i] = f;
581   }
582 
583   /**
584    * Set FreeR value flag.
585    *
586    * @param i If FreeR value is i, it is marked for cross validation.
587    */
588   public void setFreeRFlag(int i) {
589     rFreeFlag = i;
590   }
591 
592   /**
593    * setFs
594    *
595    * @param i a int.
596    * @param c a {@link ComplexNumber} object.
597    */
598   public void setFs(int i, ComplexNumber c) {
599     fs[i][0] = c.re();
600     fs[i][1] = c.im();
601   }
602 
603   /**
604    * Set amplitude sigma (sigF).
605    *
606    * @param i reflection to set
607    * @param sigF value of sigF desired
608    */
609   public void setSigF(int i, double sigF) {
610     fSigF[i][1] = sigF;
611   }
612 
613   /**
614    * setCrystalReciprocalSpace_fc
615    *
616    * @param crystalReciprocalSpace a {@link ffx.xray.CrystalReciprocalSpace} object.
617    */
618   void setCrystalReciprocalSpaceFc(CrystalReciprocalSpace crystalReciprocalSpace) {
619     this.crystalReciprocalSpaceFc = crystalReciprocalSpace;
620   }
621 
622   /**
623    * setCrystalReciprocalSpace_fs
624    *
625    * @param crystalReciprocalSpace a {@link ffx.xray.CrystalReciprocalSpace} object.
626    */
627   void setCrystalReciprocalSpaceFs(CrystalReciprocalSpace crystalReciprocalSpace) {
628     this.crystalReciprocalSpaceFs = crystalReciprocalSpace;
629   }
630 
631   /**
632    * isFreeR
633    *
634    * @param i a int.
635    * @return a boolean.
636    */
637   boolean isFreeR(int i) {
638     return (freeR[i] == rFreeFlag);
639   }
640 
641   /**
642    * get the complex number for a Fc reflection
643    *
644    * @param i reflection to get
645    * @return newly allocated {@link ComplexNumber}
646    */
647   ComplexNumber getFc(int i) {
648     return new ComplexNumber(fc[i][0], fc[i][1]);
649   }
650 
651   /**
652    * get the complex number for a Fc reflection
653    *
654    * @param i reflection to get
655    * @param c {@link ComplexNumber} to fill
656    */
657   void getFcIP(int i, ComplexNumber c) {
658     c.re(fc[i][0]);
659     c.im(fc[i][1]);
660   }
661 
662   /**
663    * getFsIP
664    *
665    * @param i a int.
666    * @param c a {@link ComplexNumber} object.
667    */
668   void getFsIP(int i, ComplexNumber c) {
669     c.re(fs[i][0]);
670     c.im(fs[i][1]);
671   }
672 
673   /**
674    * getFcTotIP
675    *
676    * @param i a int.
677    * @param c a {@link ComplexNumber} object.
678    */
679   void getFcTotIP(int i, ComplexNumber c) {
680     c.re(fcTot[i][0]);
681     c.im(fcTot[i][1]);
682   }
683 }