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 ffx.algorithms.AlgorithmListener;
41  import ffx.algorithms.Terminatable;
42  import ffx.numerics.optimization.LBFGS;
43  import ffx.numerics.optimization.LineSearch.LineSearchResult;
44  import ffx.numerics.optimization.OptimizationListener;
45  import ffx.potential.MolecularAssembly;
46  import ffx.potential.bonded.Atom;
47  import ffx.potential.bonded.Molecule;
48  import ffx.potential.bonded.Residue;
49  import ffx.realspace.RealSpaceData;
50  
51  import java.util.List;
52  import java.util.logging.Level;
53  import java.util.logging.Logger;
54  
55  import static java.lang.String.format;
56  
57  /**
58   * Refinement minimization class using {@link OptimizationListener} interface, constructs a {@link
59   * ffx.xray.RefinementEnergy} object for this purpose
60   *
61   * @author Timothy D. Fenn
62   * @since 1.0
63   */
64  public class RefinementMinimize implements OptimizationListener, Terminatable {
65  
66    private static final Logger logger = Logger.getLogger(RefinementMinimize.class.getName());
67    private static final double toSeconds = 1.0e-9;
68    public final RefinementEnergy refinementEnergy;
69    private final AlgorithmListener listener;
70    private final DataContainer dataContainer;
71    private final Atom[] activeAtomArray;
72    private final int nXYZ;
73    private final int nB;
74    private final int nOcc;
75    private final int n;
76    private final double[] x;
77    private final double[] grad;
78    private final double[] scaling;
79    private final RefinementMode refinementMode;
80    private boolean done = false;
81    private boolean terminate = false;
82    private long time;
83    private double grms;
84    private int nSteps;
85    private double eps = 0.1;
86  
87    /**
88     * constructor for refinement, assumes coordinates and B factor optimization
89     *
90     * @param data input {@link ffx.xray.DataContainer} that will be used as the model, must contain a
91     *     {@link ffx.xray.RefinementModel}
92     */
93    public RefinementMinimize(DataContainer data) {
94      this(data, RefinementMode.COORDINATES_AND_BFACTORS, null);
95    }
96  
97    /**
98     * constructor for refinement
99     *
100    * @param data input {@link ffx.xray.DataContainer} that will be used as the model, must contain a
101    *     {@link ffx.xray.RefinementModel} and either {@link ffx.xray.DiffractionData} or {@link
102    *     ffx.realspace.RealSpaceData}
103    * @param refinementmode {@link ffx.xray.RefinementMinimize.RefinementMode} for refinement
104    */
105   public RefinementMinimize(DataContainer data, RefinementMode refinementmode) {
106     this(data, refinementmode, null);
107   }
108 
109   /**
110    * constructor for refinement
111    *
112    * @param data input {@link ffx.xray.DataContainer} that will be used as the model, must contain a
113    *     {@link ffx.xray.RefinementModel} and either {@link ffx.xray.DiffractionData} or {@link
114    *     ffx.realspace.RealSpaceData}
115    * @param refinementMode {@link ffx.xray.RefinementMinimize.RefinementMode} for refinement
116    * @param listener {@link ffx.algorithms.AlgorithmListener} a listener for updates
117    */
118   public RefinementMinimize(
119       DataContainer data, RefinementMode refinementMode, AlgorithmListener listener) {
120     dataContainer = data;
121     this.listener = listener;
122     this.refinementMode = refinementMode;
123     refinementEnergy = new RefinementEnergy(data, refinementMode, null);
124     nXYZ = refinementEnergy.nXYZ;
125     nB = refinementEnergy.nBFactor;
126     nOcc = refinementEnergy.nOccupancy;
127     n = refinementEnergy.getNumberOfVariables();
128 
129     Atom[] atomArray = data.getAtomArray();
130     RefinementModel refinementModel = data.getRefinementModel();
131     int nAtoms = atomArray.length;
132 
133     // Fill an active atom array.
134     int count = 0;
135     for (Atom a : atomArray) {
136       if (a.isActive()) {
137         count++;
138       }
139     }
140     int nActive = count;
141     activeAtomArray = new Atom[count];
142     count = 0;
143     for (Atom a : atomArray) {
144       if (a.isActive()) {
145         activeAtomArray[count++] = a;
146       }
147     }
148 
149     x = new double[n];
150     grad = new double[n];
151     scaling = new double[n];
152 
153     refinementEnergy.getCoordinates(x);
154     refinementEnergy.setScaling(scaling);
155 
156     double xyzscale = 1.0;
157     double bisoscale = 1.0;
158     double anisouscale = 50.0;
159     double occscale = 15.0;
160 
161     if (refinementMode == RefinementMode.COORDINATES_AND_BFACTORS
162         || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES) {
163       bisoscale = 0.4;
164       anisouscale = 40.0;
165     }
166 
167     if (refinementMode == RefinementMode.COORDINATES_AND_OCCUPANCIES
168         || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES) {
169       occscale = 10.0;
170     }
171 
172     // set up scaling
173     if (refinementMode == RefinementMode.COORDINATES
174         || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS
175         || refinementMode == RefinementMode.COORDINATES_AND_OCCUPANCIES
176         || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES) {
177       for (int i = 0; i < n; i++) {
178         scaling[i] = xyzscale;
179       }
180     }
181 
182     if (refinementMode == RefinementMode.BFACTORS
183         || refinementMode == RefinementMode.BFACTORS_AND_OCCUPANCIES
184         || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS
185         || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES) {
186       int i = nXYZ;
187       int resnum = -1;
188       if (data instanceof DiffractionData) {
189         DiffractionData diffractionData = (DiffractionData) data;
190         int nres = diffractionData.getnResidueBFactor() + 1;
191         for (Atom a : activeAtomArray) {
192           // Ignore hydrogens or inactive atoms.
193           if (a.getAtomicNumber() == 1) {
194             continue;
195           }
196           if (a.getAnisou(null) != null) {
197             for (int j = 0; j < 6; j++) {
198               scaling[i + j] = anisouscale;
199             }
200             i += 6;
201           } else if (diffractionData.isResidueBFactor()) {
202             if (resnum != a.getResidueNumber()) {
203               if (nres >= diffractionData.getnResidueBFactor()) {
204                 if (resnum > -1 && i < nXYZ + nB - 1) {
205                   i++;
206                 }
207                 if (i < nXYZ + nB) {
208                   scaling[i] = bisoscale;
209                 }
210                 nres = 1;
211               } else {
212                 nres++;
213               }
214               resnum = a.getResidueNumber();
215             }
216           } else {
217             scaling[i] = bisoscale;
218             i++;
219           }
220         }
221       } else {
222         logger.severe(" B refinement not supported for this data type!");
223       }
224     }
225 
226     if (refinementMode == RefinementMode.OCCUPANCIES
227         || refinementMode == RefinementMode.BFACTORS_AND_OCCUPANCIES
228         || refinementMode == RefinementMode.COORDINATES_AND_OCCUPANCIES
229         || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES) {
230       if (nAtoms != nActive) {
231         logger.severe(" Occupancy refinement is not supported for inactive atoms.");
232       }
233       if (data instanceof DiffractionData) {
234 
235         int i = nXYZ + nB;
236         for (List<Residue> list : refinementModel.getAltResidues()) {
237           for (int j = 0; j < list.size(); j++) {
238             scaling[i] = occscale;
239             i++;
240           }
241         }
242         for (List<Molecule> list : refinementModel.getAltMolecules()) {
243           for (int j = 0; j < list.size(); j++) {
244             scaling[i] = occscale;
245             i++;
246           }
247         }
248       } else {
249         logger.severe(" Occupancy refinement not supported for this data type!");
250       }
251     }
252   }
253 
254   /**
255    * Parse a string into a refinement mode.
256    *
257    * @param str Refinement mode string.
258    * @return An instance of RefinementMode.
259    */
260   public static RefinementMode parseMode(String str) {
261     try {
262       return RefinementMode.valueOf(str.toUpperCase());
263     } catch (Exception e) {
264       logger.info(
265           format(
266               " Could not parse %s as a refinement mode; defaulting to coordinates.", str));
267       return RefinementMode.COORDINATES;
268     }
269   }
270 
271   /**
272    * Getter for the field <code>eps</code>.
273    *
274    * @return a {@link java.lang.Double} object.
275    */
276   public Double getEps() {
277     boolean hasaniso = false;
278 
279     for (Atom a : activeAtomArray) {
280       // ignore hydrogens!!!
281       if (a.getAtomicNumber() == 1) {
282         continue;
283       }
284       if (a.getAnisou(null) != null) {
285         hasaniso = true;
286         break;
287       }
288     }
289 
290     switch (refinementMode) {
291       case COORDINATES:
292         eps = 0.4;
293         break;
294       case BFACTORS, BFACTORS_AND_OCCUPANCIES:
295         if (hasaniso) {
296           eps = 20.0;
297         } else {
298           eps = 0.01;
299         }
300         break;
301       case COORDINATES_AND_BFACTORS, COORDINATES_AND_BFACTORS_AND_OCCUPANCIES:
302         if (hasaniso) {
303           eps = 20.0;
304         } else {
305           eps = 0.2;
306         }
307         break;
308       case OCCUPANCIES:
309         eps = 0.1;
310         break;
311       case COORDINATES_AND_OCCUPANCIES:
312         eps = 0.2;
313         break;
314     }
315 
316     return eps;
317   }
318 
319   /**
320    * get the number of B factor parameters being fit
321    *
322    * @return the number of B factor parameters
323    */
324   public int getNB() {
325     return nB;
326   }
327 
328   /**
329    * get the number of occupancy parameters being fit
330    *
331    * @return the number of occupancy parameters
332    */
333   public int getNOcc() {
334     return nOcc;
335   }
336 
337   /**
338    * get the number of xyz parameters being fit
339    *
340    * @return the number of xyz parameters
341    */
342   public int getNXYZ() {
343     return nXYZ;
344   }
345 
346   /**
347    * minimize assuming an eps of 1.0 and Integer.MAX_VALUE cycles
348    *
349    * @return {@link ffx.xray.RefinementEnergy} result
350    */
351   public RefinementEnergy minimize() {
352     return minimize(1.0);
353   }
354 
355   /**
356    * minimize assuming Integer.MAX_VALUE cycles
357    *
358    * @param eps input gradient rms desired
359    * @return {@link ffx.xray.RefinementEnergy} result
360    */
361   public RefinementEnergy minimize(double eps) {
362     return minimize(7, eps, Integer.MAX_VALUE - 2);
363   }
364 
365   /**
366    * minimize assuming an eps of 1.0 and limited cycles
367    *
368    * @param maxiter maximum iterations allowed
369    * @return {@link ffx.xray.RefinementEnergy} result
370    */
371   public RefinementEnergy minimize(int maxiter) {
372     return minimize(7, 1.0, maxiter);
373   }
374 
375   /**
376    * minimize with input eps and cycles
377    *
378    * @param eps input gradient rms desired
379    * @param maxiter maximum iterations allowed
380    * @return {@link ffx.xray.RefinementEnergy} result
381    */
382   public RefinementEnergy minimize(double eps, int maxiter) {
383     return minimize(7, eps, maxiter);
384   }
385 
386   /**
387    * minimize with input cycles for matrix conditioning, eps and cycles
388    *
389    * @param m number of cycles of matrix updates
390    * @param eps input gradient rms desired
391    * @param maxiter maximum iterations allowed
392    * @return {@link ffx.xray.RefinementEnergy} result
393    */
394   public RefinementEnergy minimize(int m, double eps, int maxiter) {
395     String typestring;
396     if (dataContainer instanceof DiffractionData) {
397       typestring = "X-ray";
398     } else if (dataContainer instanceof RealSpaceData) {
399       typestring = "Real Space";
400     } else {
401       typestring = "null";
402     }
403 
404     logger.info(" Beginning " + typestring + " Refinement");
405     switch (refinementMode) {
406       case COORDINATES:
407         logger.info(" Mode: Coordinates");
408         break;
409       case BFACTORS:
410         logger.info(" Mode: B-Factors");
411         break;
412       case COORDINATES_AND_BFACTORS:
413         logger.info(" Mode: Coordinates and B-Factors");
414         break;
415       case OCCUPANCIES:
416         logger.info(" Mode: Occupancies");
417         break;
418       case COORDINATES_AND_OCCUPANCIES:
419         logger.info(" Mode: Coordinates and Occupancies");
420         break;
421       case BFACTORS_AND_OCCUPANCIES:
422         logger.info(" Mode: B-Factors and Occupancies");
423         break;
424       case COORDINATES_AND_BFACTORS_AND_OCCUPANCIES:
425         logger.info(" Mode: Coordinates, B-Factors and Occupancies");
426         break;
427     }
428     logger.info(" Number of Parameters: " + n);
429 
430     refinementEnergy.getCoordinates(x);
431 
432     // Scale coordinates.
433     for (int i = 0; i < n; i++) {
434       x[i] *= scaling[i];
435     }
436 
437     long mtime = -System.nanoTime();
438     time = -System.nanoTime();
439     done = false;
440     int status;
441     double e = refinementEnergy.energyAndGradient(x, grad);
442     status = LBFGS.minimize(n, m, x, e, grad, eps, maxiter, refinementEnergy, this);
443     done = true;
444     switch (status) {
445       case 0:
446         logger.info(format("\n Optimization achieved convergence criteria: %8.5f", grms));
447         break;
448       case 1:
449         logger.info(format("\n Optimization terminated at step %d.", nSteps));
450         break;
451       default:
452         logger.warning("\n Optimization failed.");
453     }
454 
455     if (logger.isLoggable(Level.INFO)) {
456       StringBuilder sb = new StringBuilder();
457       mtime += System.nanoTime();
458       sb.append(format(" Optimization time: %g (sec)", mtime * toSeconds));
459       logger.info(sb.toString());
460     }
461 
462     return refinementEnergy;
463   }
464 
465   /** {@inheritDoc} */
466   @Override
467   public boolean optimizationUpdate(
468       int iter,
469       int nBFGS,
470       int nfun,
471       double grms,
472       double xrms,
473       double f,
474       double df,
475       double angle,
476       LineSearchResult info) {
477     long currentTime = System.nanoTime();
478     Double seconds = (currentTime - time) * 1.0e-9;
479     time = currentTime;
480     this.grms = grms;
481     this.nSteps = iter;
482 
483     // Update display.
484     if (listener != null) {
485       MolecularAssembly[] molecularAssembly = dataContainer.getMolecularAssemblies();
486       for (MolecularAssembly ma : molecularAssembly) {
487         listener.algorithmUpdate(ma);
488       }
489     }
490 
491     if (iter == 0) {
492       if (nBFGS > 0) {
493         logger.info("\n Limited Memory BFGS Quasi-Newton Optimization: \n");
494       } else {
495         logger.info("\n Steepest Decent Optimization: \n");
496       }
497       logger.info(
498           " Cycle       Energy      G RMS    Delta E   Delta X    Angle  Evals     Time      "
499               + dataContainer.printOptimizationHeader());
500     }
501     if (info == null) {
502       logger.info(format("%6d %12.3f %10.3f", iter, f, grms));
503     } else if (info == LineSearchResult.Success) {
504       StringBuilder sb = new StringBuilder();
505       sb.append(format("%6d %12.3f %10.3f %10.3f %9.4f %8.2f %6d %8.3f ",
506               iter, f, grms, df, xrms, angle, nfun, seconds));
507       sb.append(dataContainer.printOptimizationUpdate());
508       logger.info(sb.toString());
509     } else {
510       logger.info(
511           format("%6d %12.3f %10.3f %10.3f %9.4f %8.2f %6d %8s",
512               iter, f, grms, df, xrms, angle, nfun, info));
513     }
514     if (terminate) {
515       logger.info(" The optimization recieved a termination request.");
516       // Tell the L-BFGS optimizer to terminate.
517       return false;
518     }
519     return true;
520   }
521 
522   /** {@inheritDoc} */
523   @Override
524   public void terminate() {
525     terminate = true;
526     while (!done) {
527       synchronized (this) {
528         try {
529           wait(1);
530         } catch (Exception e) {
531           logger.log(Level.WARNING, " Exception terminating minimization.\n", e);
532         }
533       }
534     }
535   }
536 
537   /** Different refinement mode selection types. */
538   public enum RefinementMode {
539 
540     /** refine coordinates only */
541     COORDINATES,
542     /** refine B factors only (if anisotropic, refined as such) */
543     BFACTORS,
544     /** refine coordinates and B factors (if anisotropic, refined as such) */
545     COORDINATES_AND_BFACTORS,
546     /** refine occupancies only (alternate conformers are constrained) */
547     OCCUPANCIES,
548     /** refine B factors and occupancies */
549     BFACTORS_AND_OCCUPANCIES,
550     /** refine coordinates and occupancies */
551     COORDINATES_AND_OCCUPANCIES,
552     /** refine all */
553     COORDINATES_AND_BFACTORS_AND_OCCUPANCIES
554   }
555 }