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.algorithms.optimize;
39  
40  import ffx.algorithms.AlgorithmListener;
41  import ffx.numerics.Potential;
42  import ffx.numerics.optimization.LineSearch;
43  import ffx.potential.ForceFieldEnergy;
44  import ffx.potential.openmm.OpenMMEnergy;
45  import ffx.potential.MolecularAssembly;
46  import ffx.potential.openmm.OpenMMContext;
47  import ffx.potential.openmm.OpenMMState;
48  import ffx.potential.bonded.Atom;
49  
50  import java.util.logging.Logger;
51  
52  import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_DataType.OpenMM_State_Energy;
53  import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_DataType.OpenMM_State_Forces;
54  import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_DataType.OpenMM_State_Positions;
55  import static java.lang.Double.isInfinite;
56  import static java.lang.Double.isNaN;
57  import static java.lang.String.format;
58  import static org.apache.commons.math3.util.FastMath.sqrt;
59  
60  /**
61   * OpenMM accelerated L-BFGS minimization.
62   *
63   * @author Michael J. Schnieders
64   * @since 1.0
65   */
66  public class MinimizeOpenMM extends Minimize {
67  
68    private static final Logger logger = Logger.getLogger(MinimizeOpenMM.class.getName());
69  
70    public MinimizeOpenMM(MolecularAssembly molecularAssembly) {
71      super(molecularAssembly, molecularAssembly.getPotentialEnergy(), null);
72    }
73  
74    public MinimizeOpenMM(MolecularAssembly molecularAssembly,
75                          OpenMMEnergy openMMEnergy) {
76      super(molecularAssembly, openMMEnergy, null);
77    }
78  
79    public MinimizeOpenMM(MolecularAssembly molecularAssembly,
80                          OpenMMEnergy openMMEnergy, AlgorithmListener algorithmListener) {
81      super(molecularAssembly, openMMEnergy, algorithmListener);
82    }
83  
84    /**
85     * Note the OpenMM L-BFGS minimizer does not accept the parameter "m" for the number of previous
86     * steps used to estimate the Hessian.
87     *
88     * @param m             The number of previous steps used to estimate the Hessian (ignored).
89     * @param eps           The convergence criteria.
90     * @param maxIterations The maximum number of iterations.
91     * @return The potential.
92     */
93    @Override
94    public Potential minimize(int m, double eps, int maxIterations) {
95      return minimize(eps, maxIterations);
96    }
97  
98    /**
99     * minimize
100    *
101    * @param eps           The convergence criteria.
102    * @param maxIterations The maximum number of iterations.
103    * @return a {@link ffx.numerics.Potential} object.
104    */
105   @Override
106   public Potential minimize(double eps, int maxIterations) {
107 
108     ForceFieldEnergy forceFieldEnergy = molecularAssembly.getPotentialEnergy();
109 
110     if (forceFieldEnergy instanceof OpenMMEnergy openMMEnergy) {
111       time = -System.nanoTime();
112 
113       // Respect the use flag, and lambda state.
114       Atom[] atoms = molecularAssembly.getAtomArray();
115       openMMEnergy.updateParameters(atoms);
116 
117       // Respect (in)active atoms.
118       openMMEnergy.setActiveAtoms();
119 
120       // Get the coordinates to start from.
121       openMMEnergy.getCoordinates(x);
122 
123       // Calculate the starting energy before optimization.
124       double e = openMMEnergy.energy(x);
125       logger.info(format("\n Initial energy:                 %12.6f (kcal/mol)", e));
126 
127       // Run the minimization in the current OpenMM Context.
128       OpenMMContext openMMContext = openMMEnergy.getContext();
129       openMMContext.optimize(eps, maxIterations);
130 
131       // Get the minimized coordinates, forces and potential energy back from OpenMM.
132       int mask = OpenMM_State_Energy | OpenMM_State_Positions | OpenMM_State_Forces;
133       OpenMMState openMMState = openMMContext.getOpenMMState(mask);
134       energy = openMMState.potentialEnergy;
135       openMMState.getPositions(x);
136       openMMState.getGradient(grad);
137       openMMState.destroy();
138 
139       // Compute the RMS gradient.
140       int index = 0;
141       double grad2 = 0;
142       for (Atom atom : atoms) {
143         if (atom.isActive()) {
144           double fx = grad[index++];
145           double fy = grad[index++];
146           double fz = grad[index++];
147           grad2 += fx * fx + fy * fy + fz * fz;
148         }
149       }
150       rmsGradient = sqrt(grad2 / n);
151 
152       double[] ffxGrad = new double[n];
153       openMMEnergy.getCoordinates(x);
154       double ffxEnergy = openMMEnergy.energyAndGradientFFX(x, ffxGrad);
155       double grmsFFX = 0.0;
156       for (int i = 0; i < n; i++) {
157         double gi = ffxGrad[i];
158         if (isNaN(gi) || isInfinite(gi)) {
159           String message = format(" The gradient of variable %d is %8.3f.", i, gi);
160           logger.warning(message);
161         }
162         grmsFFX += gi * gi;
163       }
164       grmsFFX = sqrt(grmsFFX / n);
165 
166       time += System.nanoTime();
167       logger.info(format(" Final energy for OpenMM         %12.6f vs. FFX %12.6f in %8.3f (sec).", energy, ffxEnergy, time * 1.0e-9));
168       logger.info(format(" Convergence criteria for OpenMM %12.6f vs. FFX %12.6f (kcal/mol/A).", rmsGradient, grmsFFX));
169     }
170 
171     if (algorithmListener != null) {
172       algorithmListener.algorithmUpdate(molecularAssembly);
173     }
174 
175     return forceFieldEnergy;
176   }
177 
178   /**
179    * MinimizeOpenMM does not support the OptimizationListener interface.
180    *
181    * @since 1.0
182    */
183   @Override
184   public boolean optimizationUpdate(int iteration, int nBFGS, int functionEvaluations,
185                                     double rmsGradient, double rmsCoordinateChange, double energy, double energyChange,
186                                     double angle, LineSearch.LineSearchResult lineSearchResult) {
187     logger.warning(" MinimizeOpenMM does not support updates at each optimization step.");
188     return false;
189   }
190 }