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;
39  
40  import static java.lang.Double.isNaN;
41  import static java.lang.String.format;
42  import static java.lang.System.arraycopy;
43  import static java.util.Arrays.fill;
44  import static org.apache.commons.math3.util.FastMath.log;
45  
46  import edu.rit.pj.ParallelTeam;
47  import ffx.algorithms.Terminatable;
48  import ffx.crystal.Crystal;
49  import ffx.crystal.HKL;
50  import ffx.crystal.ReflectionList;
51  import ffx.numerics.optimization.LBFGS;
52  import ffx.numerics.optimization.LineSearch.LineSearchResult;
53  import ffx.numerics.optimization.OptimizationListener;
54  import ffx.xray.CrystalReciprocalSpace.SolventModel;
55  
56  import javax.annotation.Nullable;
57  import java.util.logging.Level;
58  import java.util.logging.Logger;
59  
60  /**
61   * ScaleBulkMinimize class.
62   *
63   * @author Timothy D. Fenn
64   * @since 1.0
65   */
66  public class ScaleBulkMinimize implements OptimizationListener, Terminatable {
67  
68    private static final Logger logger = Logger.getLogger(ScaleBulkMinimize.class.getName());
69    private static final double toSeconds = 1.0e-9;
70    private final ReflectionList reflectionlist;
71    private final DiffractionRefinementData refinementData;
72    private final Crystal crystal;
73    private final CrystalReciprocalSpace crystalReciprocalSpace;
74    private final ScaleBulkEnergy bulkSolventEnergy;
75    private final int solventN;
76    private final int n;
77    private final double[] x;
78    private final double[] grad;
79    private final double[] scaling;
80    private boolean done = false;
81    private boolean terminate = false;
82    private long time;
83    private double grms;
84    private int nSteps;
85  
86    /**
87     * Constructor for ScaleBulkMinimize.
88     *
89     * @param reflectionList a {@link ffx.crystal.ReflectionList} object.
90     * @param refinementData a {@link ffx.xray.DiffractionRefinementData} object.
91     * @param crystalReciprocalSpace a {@link ffx.xray.CrystalReciprocalSpace} object.
92     * @param parallelTeam the ParallelTeam to execute the ScaleBulkMinimize.
93     */
94    public ScaleBulkMinimize(
95        ReflectionList reflectionList,
96        DiffractionRefinementData refinementData,
97        CrystalReciprocalSpace crystalReciprocalSpace,
98        ParallelTeam parallelTeam) {
99      this.reflectionlist = reflectionList;
100     this.refinementData = refinementData;
101     this.crystal = reflectionList.crystal;
102     this.crystalReciprocalSpace = crystalReciprocalSpace;
103 
104     if (crystalReciprocalSpace.solventModel == SolventModel.NONE) {
105       solventN = 1;
106     } else {
107       solventN = 3;
108     }
109     n = solventN + refinementData.nScale;
110     bulkSolventEnergy = new ScaleBulkEnergy(reflectionList, refinementData, n, parallelTeam);
111 
112     x = new double[n];
113     grad = new double[n];
114     scaling = new double[n];
115     fill(scaling, 1.0);
116 
117     x[0] = refinementData.modelScaleK;
118     if (solventN > 1) {
119       x[1] = refinementData.bulkSolventK;
120       x[2] = refinementData.bulkSolventUeq;
121     }
122     for (int i = 0; i < 6; i++) {
123       if (crystal.scaleB[i] >= 0) {
124         x[solventN + crystal.scaleB[i]] = refinementData.modelAnisoB[i];
125       }
126     }
127 
128     setInitialScale();
129   }
130 
131   /**
132    * getCoordinates.
133    *
134    * @param x the array to populate with parameters or null to create a new array.
135    * @return an array containing the parameters.
136    */
137   public double[] getCoordinates(@Nullable double[] x) {
138     if (x == null) {
139       x = new double[this.x.length];
140     }
141     arraycopy(this.x, 0, x, 0, this.x.length);
142     return x;
143   }
144 
145   /**
146    * getNumberOfVariables.
147    *
148    * @return a int.
149    */
150   public int getNumberOfVariables() {
151     return x.length;
152   }
153 
154   /** ksbsGridOptimize */
155   public void ksbsGridOptimize() {
156     if (solventN < 3) {
157       return;
158     }
159 
160     bulkSolventEnergy.setScaling(scaling);
161 
162     double min = Double.POSITIVE_INFINITY;
163     double bulkSolventK = refinementData.bulkSolventK;
164     double bulkSolventUeq = refinementData.bulkSolventUeq;
165     double kmin = 0.05;
166     double kmax = 0.9;
167     double kstep = 0.05;
168     double bmin = 10.0;
169     double bmax = 200.0;
170     double bstep = 10.0;
171     for (double i = kmin; i <= kmax; i += kstep) {
172       for (double j = bmin; j <= bmax; j += bstep) {
173         x[1] = i;
174         x[2] = j;
175         double sum = bulkSolventEnergy.energyAndGradient(x, grad);
176         logger.info(" ks: " + i + " bs: " + j + " sum: " + sum);
177         if (sum < min) {
178           min = sum;
179           bulkSolventK = i;
180           bulkSolventUeq = j;
181         }
182       }
183     }
184     logger.info(" minks: " + bulkSolventK + " minbs: " + bulkSolventUeq + " min: " + min);
185     refinementData.bulkSolventK = bulkSolventK;
186     refinementData.bulkSolventUeq = bulkSolventUeq;
187     bulkSolventEnergy.setScaling(null);
188   }
189 
190   /**
191    * minimize
192    *
193    * @return a {@link ffx.xray.ScaleBulkEnergy} object.
194    */
195   public ScaleBulkEnergy minimize() {
196     return minimize(0.5);
197   }
198 
199   /**
200    * minimize
201    *
202    * @param eps a double.
203    * @return a {@link ffx.xray.ScaleBulkEnergy} object.
204    */
205   public ScaleBulkEnergy minimize(double eps) {
206     return minimize(5, eps);
207   }
208 
209   /**
210    * minimize
211    *
212    * @param m a int.
213    * @param eps a double.
214    * @return a {@link ffx.xray.ScaleBulkEnergy} object.
215    */
216   public ScaleBulkEnergy minimize(int m, double eps) {
217     bulkSolventEnergy.setScaling(scaling);
218     double e = bulkSolventEnergy.energyAndGradient(x, grad);
219 
220     long mtime = -System.nanoTime();
221     time = -System.nanoTime();
222     done = false;
223     int status = LBFGS.minimize(n, m, x, e, grad, eps, bulkSolventEnergy, this);
224     done = true;
225     switch (status) {
226       case 0:
227         logger.info(format("\n Optimization achieved convergence criteria: %10.5f\n", grms));
228         break;
229       case 1:
230         logger.info(format("\n Optimization terminated at step %d.\n", nSteps));
231         break;
232       default:
233         logger.warning("\n Optimization failed.\n");
234     }
235     refinementData.modelScaleK = x[0] / scaling[0];
236     if (solventN > 1) {
237       refinementData.bulkSolventK = x[1] / scaling[1];
238       refinementData.bulkSolventUeq = x[2] / scaling[2];
239       if (crystalReciprocalSpace != null) {
240         refinementData.solventA = crystalReciprocalSpace.solventA;
241         refinementData.solventB = crystalReciprocalSpace.solventB;
242       }
243     }
244     for (int i = 0; i < 6; i++) {
245       int offset = crystal.scaleB[i];
246       if (offset >= 0) {
247         int index = solventN + offset;
248         refinementData.modelAnisoB[i] = x[index] / scaling[index];
249       }
250     }
251 
252     if (logger.isLoggable(Level.INFO)) {
253       mtime += System.nanoTime();
254       logger.info(format(" Optimization time: %8.3f (sec)\n", mtime * toSeconds));
255     }
256 
257     bulkSolventEnergy.setScaling(null);
258 
259     return bulkSolventEnergy;
260   }
261 
262   /** {@inheritDoc} */
263   @Override
264   public boolean optimizationUpdate(
265       int iter,
266       int nBFGS,
267       int nfun,
268       double grms,
269       double xrms,
270       double f,
271       double df,
272       double angle,
273       LineSearchResult info) {
274     long currentTime = System.nanoTime();
275     Double seconds = (currentTime - time) * 1.0e-9;
276     time = currentTime;
277     this.grms = grms;
278     this.nSteps = iter;
279 
280     if (iter == 0) {
281       if (nBFGS > 0) {
282         logger.info("\n Limited Memory BFGS Quasi-Newton Optimization of Fc to Fo Scaling Parameters\n");
283       } else {
284         logger.info("\n Steepest Decent Optimization of Fc to Fo Scaling Parameters\n");
285       }
286       logger.info(" Cycle       Energy      G RMS    Delta E   Delta X    Angle  Evals     Time");
287     }
288     if (info == null) {
289       logger.info(format("%6d %12.5f %10.6f", iter, f, grms));
290     } else {
291       if (info == LineSearchResult.Success) {
292         logger.info(format("%6d %12.5f %10.6f %10.6f %9.5f %8.2f %6d %8.3f",
293                 iter, f, grms, df, xrms, angle, nfun, seconds));
294       } else {
295         logger.info(format("%6d %12.5f %10.6f %10.6f %9.5f %8.2f %6d %8s",
296                 iter, f, grms, df, xrms, angle, nfun, info));
297       }
298     }
299     if (terminate) {
300       logger.info(" The optimization recieved a termination request.");
301       // Tell the L-BFGS optimizer to terminate.
302       return false;
303     }
304     return true;
305   }
306 
307   /** {@inheritDoc} */
308   @Override
309   public void terminate() {
310     terminate = true;
311     while (!done) {
312       synchronized (this) {
313         try {
314           wait(1);
315         } catch (Exception e) {
316           logger.log(Level.WARNING, "Exception terminating minimization.\n", e);
317         }
318       }
319     }
320   }
321 
322   /**
323    * getScaleBulkEnergy.
324    *
325    * @return a {@link ffx.xray.ScaleBulkEnergy} object.
326    */
327   ScaleBulkEnergy getScaleBulkEnergy() {
328     return bulkSolventEnergy;
329   }
330 
331   private void setInitialScale() {
332     double[][] fc = refinementData.fc;
333     double[][] fSigf = refinementData.fSigF;
334 
335     bulkSolventEnergy.setScaling(scaling);
336     double e = bulkSolventEnergy.energyAndGradient(x, grad);
337     bulkSolventEnergy.setScaling(null);
338 
339     double sumfofc = 0.0;
340     double sumfc = 0.0;
341     for (HKL ih : reflectionlist.hklList) {
342       int i = ih.getIndex();
343       if (isNaN(fc[i][0]) || isNaN(fSigf[i][0]) || fSigf[i][1] <= 0.0) {
344         continue;
345       }
346 
347       double fcTotF = refinementData.fcTotF(i);
348       sumfofc += fSigf[i][0] * fcTotF;
349       sumfc += fcTotF * fcTotF;
350     }
351 
352     x[0] = log(4.0 * sumfofc / sumfc);
353   }
354 
355   /** GridOptimize */
356   void gridOptimize() {
357     if (crystalReciprocalSpace == null) {
358       return;
359     }
360 
361     bulkSolventEnergy.setScaling(scaling);
362 
363     double min = Double.POSITIVE_INFINITY;
364     double solventA = crystalReciprocalSpace.solventA;
365     double solventB = crystalReciprocalSpace.solventB;
366     double amin = solventA - 1.0;
367     double amax = (solventA + 1.0) / 0.9999;
368     double astep = 0.25;
369     double bmin = solventB - 0.2;
370     double bmax = (solventB + 0.2) / 0.9999;
371     double bstep = 0.05;
372     if (crystalReciprocalSpace.solventModel == SolventModel.BINARY) {
373       amin = solventA - 0.2;
374       amax = (solventA + 0.2) / 0.9999;
375       astep = 0.05;
376     }
377 
378     logger.info(" Bulk Solvent Grid Search");
379     for (double i = amin; i <= amax; i += astep) {
380       for (double j = bmin; j <= bmax; j += bstep) {
381         crystalReciprocalSpace.setSolventAB(i, j);
382         crystalReciprocalSpace.computeDensity(refinementData.fs);
383         double sum = bulkSolventEnergy.energy(x);
384         logger.info(format(" A: %6.3f B: %6.3f Sum: %12.8f", i, j, sum));
385         if (sum < min) {
386           min = sum;
387           solventA = i;
388           solventB = j;
389         }
390       }
391     }
392 
393     logger.info(format("\n Minimum at\n A: %6.3f B: %6.3f Sum: %12.8f", solventA, solventB, min));
394     crystalReciprocalSpace.setSolventAB(solventA, solventB);
395     refinementData.solventA = solventA;
396     refinementData.solventB = solventB;
397     crystalReciprocalSpace.computeDensity(refinementData.fs);
398     bulkSolventEnergy.setScaling(null);
399   }
400 }