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.algorithms.optimize.anneal;
39  
40  import ffx.algorithms.AlgorithmListener;
41  import ffx.algorithms.Terminatable;
42  import ffx.algorithms.dynamics.MolecularDynamics;
43  import ffx.algorithms.dynamics.integrators.IntegratorEnum;
44  import ffx.algorithms.dynamics.thermostats.ThermostatEnum;
45  import ffx.numerics.Potential;
46  import ffx.potential.MolecularAssembly;
47  
48  import java.io.File;
49  import java.util.Set;
50  import java.util.logging.Level;
51  import java.util.logging.Logger;
52  
53  /**
54   * Run NVT molecular dynamics at a series of temperatures to optimize a structure.
55   *
56   * @author Michael J. Schnieders
57   * @author Jacob M. Litman
58   * @since 1.0
59   */
60  public class SimulatedAnnealing implements Runnable, Terminatable {
61  
62    private static final Logger logger = Logger.getLogger(SimulatedAnnealing.class.getName());
63    /** The MolecularDynamics instance used for Simulated Annealing. */
64    private final MolecularDynamics molecularDynamics;
65    /** Schedule for annealing. */
66    private final AnnealingSchedule schedule;
67    /** Number of MD steps per annealing window. */
68    private final long mdSteps;
69    /** Integration time step. */
70    private final double timeStep;
71    /** Whether to reinitialize velocities at the start of each time step. */
72    private final boolean reInitVelocity;
73    /** Number of MD steps per OpenMM cycle (assuming OpenMM is used!). */
74    private final int trajSteps = 1;
75    /** Interval to print updates to the screen. */
76    private double printInterval = 0.01;
77    /** Flag to indicate the algorithm is done. */
78    private boolean done = true;
79    /** Flag to indicate the UI has requested the algorithm terminate. */
80    private boolean terminate;
81  
82    private double saveFrequency = 0.1;
83    /** Restart file. */
84    private File dynFile;
85  
86    /**
87     * Constructor for SimulatedAnnealing.
88     *
89     * @param molecularAssembly The Molecular Assembly to operate on.
90     * @param potentialEnergy The potential to anneal against.
91     * @param algorithmListener The algorithm listener is a callback to UI.
92     * @param requestedThermostat The requested thermostat.
93     * @param requestedIntegrator The requested integrator.
94     * @param annealingSchedule Schedule of temperatures to simulate at.
95     * @param mdSteps Steps per SA window.
96     * @param timeStep Time step for MD in psec.
97     * @param reInitVelocity Request velocities to be reinitialized.
98     * @param dynFile Dynamics restart file to begin from.
99     */
100   public SimulatedAnnealing(MolecularAssembly molecularAssembly, Potential potentialEnergy,
101       AlgorithmListener algorithmListener, ThermostatEnum requestedThermostat,
102       IntegratorEnum requestedIntegrator, AnnealingSchedule annealingSchedule, long mdSteps,
103       double timeStep, boolean reInitVelocity, File dynFile) {
104 
105     molecularDynamics = MolecularDynamics.dynamicsFactory(molecularAssembly, potentialEnergy,
106         algorithmListener, requestedThermostat, requestedIntegrator);
107     this.schedule = annealingSchedule;
108     this.mdSteps = mdSteps;
109     this.timeStep = timeStep;
110     this.reInitVelocity = reInitVelocity;
111     this.dynFile = dynFile;
112   }
113 
114   /** anneal */
115   public void anneal() {
116     // Return if already running; Could happen if two threads call anneal
117     // on the same SimulatedAnnealing instance.
118     if (!done) {
119       logger.warning(" Programming error - a thread invoked anneal when it was already running.");
120       return;
121     }
122     done = false;
123     logger.info(" Beginning simulated annealing");
124     begin();
125   }
126 
127   /**
128    * getKineticEnergy.
129    *
130    * @return a double.
131    */
132   public double getKineticEnergy() {
133     return molecularDynamics.getKineticEnergy();
134   }
135 
136   /**
137    * getPotentialEnergy.
138    *
139    * @return a double.
140    */
141   public double getPotentialEnergy() {
142     return molecularDynamics.getPotentialEnergy();
143   }
144 
145   /**
146    * getTemperature.
147    *
148    * @return a double.
149    */
150   public double getTemperature() {
151     return molecularDynamics.getTemperature();
152   }
153 
154   /**
155    * getTotalEnergy.
156    *
157    * @return a double.
158    */
159   public double getTotalEnergy() {
160     return molecularDynamics.getTotalEnergy();
161   }
162 
163   /**
164    * {@inheritDoc}
165    *
166    * <p>This method should only be invoked within the SimulatedAnnealing instance.
167    */
168   @Override
169   public void run() {
170     done = false;
171     terminate = false;
172 
173     int minMdSteps = (int) (mdSteps * schedule.minWindowLength());
174     if (minMdSteps < trajSteps) {
175       logger.warning(String.format(
176           " Minimum number of MD steps per annealing cycle %d was less than steps per OpenMM MD cycle %d! Setting steps per MD cycle to %d",
177           minMdSteps, trajSteps, minMdSteps));
178       setTrajectorySteps(minMdSteps);
179     }
180 
181     int nWindows = schedule.getNumWindows();
182     boolean forceFirstReinit = (dynFile == null);
183 
184     for (int i = 0; i < nWindows; i++) {
185       double temperature = schedule.getTemperature(i);
186       int nSteps = (int) (schedule.windowLength(i) * mdSteps);
187       logger.info(
188           String.format(" Annealing window %d: %d steps at %9.4g K", (i + 1), nSteps, temperature));
189       molecularDynamics.dynamic(nSteps, timeStep, printInterval, saveFrequency, temperature,
190           (reInitVelocity || forceFirstReinit), dynFile);
191       if (dynFile == null) {
192         dynFile = molecularDynamics.getDynFile();
193       }
194       forceFirstReinit = false;
195       if (terminate) {
196         logger.info(String.format("\n Terminating at temperature %8.3f.\n", temperature));
197         break;
198       }
199     }
200     if (!terminate) {
201       logger.info(String.format(" Completed %8d annealing steps\n", nWindows));
202     }
203 
204     done = true;
205     terminate = false;
206   }
207 
208   /**
209    * Setter for the field <code>printInterval</code>.
210    *
211    * @param printInterval a double.
212    */
213   public void setPrintInterval(double printInterval) {
214     this.printInterval = printInterval;
215   }
216 
217   /**
218    * Method to set the Restart Frequency.
219    *
220    * @param restart the time between writing restart files.
221    * @throws java.lang.IllegalArgumentException If restart frequency is not a positive number
222    */
223   public void setRestartFrequency(double restart) throws IllegalArgumentException {
224     if (Double.isFinite(restart) && restart > 0) {
225       molecularDynamics.setRestartFrequency(restart);
226     } else {
227       throw new IllegalArgumentException(
228           String.format(" Restart frequency must be positive finite, was %10.4g", restart));
229     }
230   }
231 
232   /**
233    * Sets the frequency of writing to the trajectory file.
234    *
235    * @param save Frequency (psec^-1) to write out the trajectory.
236    */
237   public void setSaveFrequency(double save) {
238     this.saveFrequency = save;
239   }
240 
241   /**
242    * Sets the number of steps to use per OpenMM cycle.
243    *
244    * @param trajectorySteps Steps per OpenMM cycle.
245    */
246   public void setTrajectorySteps(int trajectorySteps) {
247     molecularDynamics.setIntervalSteps(trajectorySteps);
248   }
249 
250   /** {@inheritDoc} */
251   @Override
252   public void terminate() {
253     terminate = true;
254     while (!done) {
255       synchronized (this) {
256         try {
257           wait(1);
258         } catch (Exception e) {
259           logger.log(Level.WARNING, "Exception terminating annealing.\n", e);
260         }
261       }
262     }
263   }
264 
265   private void begin() {
266     logger.info(String.format(" Initial temperature:    %8.3f (Kelvin)", schedule.getHighTemp()));
267     logger.info(String.format(" Final temperature:      %8.3f (Kelvin)", schedule.getLowTemp()));
268     logger.info(String.format(" Annealing steps:        %8d", schedule.getNumWindows()));
269     logger.info(String.format(" MD steps/temperature:   %8d", mdSteps));
270     logger.info(String.format(" MD time step:           %8.3f (fs)", timeStep));
271 
272     Thread annealingThread = new Thread(this);
273     annealingThread.start();
274     synchronized (this) {
275       try {
276         while (annealingThread.isAlive()) {
277           wait(100);
278         }
279       } catch (Exception e) {
280         String message = "Simulated annealing interrupted.";
281         logger.log(Level.WARNING, message, e);
282       }
283     }
284   }
285 
286   /** Represents non-composite AnnealingSchedules known (i.e. not FlatEndAnnealSchedule). */
287   public enum Schedules {
288     EXP(ExpAnnealSchedule::new, "EXP", "EXPONENTIAL"), LINEAR(LinearAnnealSchedule::new, "LINEAR");
289 
290     private final ScheduleConstructor sc;
291     private final Set<String> aliases;
292 
293     Schedules(ScheduleConstructor sc, String... names) {
294       this.sc = sc;
295       aliases = Set.of(names);
296     }
297 
298     /**
299      * Attempt to parse a String to a Schedules in a case-insensitive, alias-recognizing fashion.
300      *
301      * @param name Name of a schedule.
302      * @return A Schedules enum.
303      */
304     public static Schedules parse(String name) {
305       name = name.toUpperCase();
306       for (Schedules s : values()) {
307         if (s.aliases.contains(name)) {
308           return s;
309         }
310       }
311       return valueOf(name);
312     }
313 
314     /**
315      * Creates an AnnealingSchedule corresponding to this enum and provided values.
316      *
317      * @param nWindows Number of annealing windows.
318      * @param tLow Final temperature.
319      * @param tHigh Starting temperature.
320      * @return An AnnealingSchedule.
321      */
322     public AnnealingSchedule generate(int nWindows, double tLow, double tHigh) {
323       return sc.asConstruct(nWindows, tLow, tHigh);
324     }
325   }
326 
327   /** Functional interface corresponding to constructors of non-composite AnnealingSchedules. */
328   private interface ScheduleConstructor {
329 
330     AnnealingSchedule asConstruct(int nWindows, double tLow, double tHigh);
331   }
332 }