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.thermodynamics;
39  
40  import ffx.algorithms.AlgorithmListener;
41  import ffx.algorithms.cli.DynamicsOptions;
42  import ffx.algorithms.dynamics.MDVerbosity;
43  import ffx.algorithms.dynamics.MDWriteAction;
44  import ffx.algorithms.dynamics.MolecularDynamics;
45  import ffx.algorithms.dynamics.integrators.IntegratorEnum;
46  import ffx.algorithms.dynamics.thermostats.ThermostatEnum;
47  import ffx.algorithms.mc.BoltzmannMC;
48  import ffx.algorithms.mc.LambdaMove;
49  import ffx.algorithms.mc.MDMove;
50  import ffx.algorithms.thermodynamics.OrthogonalSpaceTempering.Histogram;
51  import ffx.numerics.Potential;
52  import ffx.potential.MolecularAssembly;
53  import ffx.potential.utils.EnergyException;
54  import org.apache.commons.configuration2.CompositeConfiguration;
55  
56  import java.io.File;
57  import java.util.EnumSet;
58  import java.util.logging.Level;
59  import java.util.logging.Logger;
60  
61  import static ffx.utilities.Constants.NS2SEC;
62  import static java.lang.String.format;
63  import static java.lang.System.arraycopy;
64  import static java.lang.System.nanoTime;
65  import static org.apache.commons.math3.util.FastMath.abs;
66  
67  /**
68   * Sample a thermodynamic path using the OST method, with the time-dependent bias built up using
69   * Metropolis Monte Carlo steps.
70   *
71   * <p>The algorithm generates coordinate (X) MC moves using molecular dynamics at a fixed lambda
72   * value (i.e., using OpenMM), followed by MC lambda moves.
73   *
74   * <p>1.) At a fixed Lambda, run a defined length MD trajectory to "move" coordinates and dU/dL on
75   * an approximate potential U* (i.e., no OST Bias).
76   *
77   * <p>2.) Accept / Reject the MD move with probability exp[-Beta(dU - dU*)] where dU is the change
78   * in AMOEBA + Bias energy and dU* is the change in AMOEBA + Kinetic energy from the MD.
79   *
80   * <p>3.) Randomly change the value of Lambda.
81   *
82   * <p>4.) Accept / Reject the Lambda move using the AMOEBA + OST Bias energy.
83   *
84   * <p>5.) Add to the time dependent 2D bias using the current values of Lambda and dU/dL.
85   *
86   * @author Michael J. Schnieders
87   * @author Hernan Beranbe
88   * @author Mallory R. Tollefson
89   * @author Jacob Litman
90   * @since 1.0
91   */
92  public class MonteCarloOST extends BoltzmannMC {
93  
94    /**
95     * Logger object to print out information for this class.
96     */
97    private static final Logger logger = Logger.getLogger(MonteCarloOST.class.getName());
98  
99    /**
100    * Controls the effect of verbose by logging at FINE vs. INFO.
101    */
102   private final Level mcLogLevel;
103   /**
104    * How verbose MD should be.
105    */
106   private final MDVerbosity mdVerbosity;
107 
108   /**
109    * The MD moves are only valid if energy is conserved. For this reason, energy drift is monitored.
110    */
111   private static final double ENERGY_CONSERVATION_TOLERANCE = 10.0;
112   /**
113    * Potential object used to retrieve the coordinates for the system.
114    */
115   private final Potential potential;
116   /**
117    * OST object used to retrieve OST energy throughout the simulation.
118    */
119   private final OrthogonalSpaceTempering orthogonalSpaceTempering;
120   /**
121    * MDMove object for completing MC-OST molecular dynamics moves.
122    */
123   private final MDMove mdMove;
124   /**
125    * Deposit a bias once every N MC cycles. Defaults to 1.
126    */
127   private final int biasDepositionFrequency;
128   /**
129    * Total number of steps to take for MC-OST sampling.
130    */
131   private long totalSteps;
132   /**
133    * Number of MD steps to take per MC-OST round.
134    */
135   private final long stepsPerMove;
136   /**
137    * Lambda move object for completing MC-OST lambda moves.
138    */
139   private final LambdaMove lambdaMove;
140   /**
141    * Double that keeps track of our lambda value.
142    */
143   private double lambda = 1.0;
144   /**
145    * Boolean that tells algorithm that we are in the equilibration phase of MC-OST.
146    */
147   private boolean equilibration = false;
148   /**
149    * True if MC-OST should handle writing out files.
150    */
151   private boolean automaticWriteouts = true;
152 
153   /**
154    * Constructor for MonteCarloOST.
155    *
156    * @param potentialEnergy          a {@link ffx.numerics.Potential} object.
157    * @param orthogonalSpaceTempering a {@link OrthogonalSpaceTempering} object.
158    * @param molecularAssembly        a {@link ffx.potential.MolecularAssembly} object.
159    * @param properties               a {@link org.apache.commons.configuration2.CompositeConfiguration} object.
160    * @param listener                 a {@link ffx.algorithms.AlgorithmListener} object.
161    * @param dynamics                 CLI object containing key information.
162    * @param verbose                  Whether to be verbose.
163    * @param cycleLength              Length of an MC cycle in MD steps.
164    */
165   public MonteCarloOST(Potential potentialEnergy, OrthogonalSpaceTempering orthogonalSpaceTempering,
166                        MolecularAssembly molecularAssembly, CompositeConfiguration properties,
167                        AlgorithmListener listener, DynamicsOptions dynamics,
168                        boolean verbose, int cycleLength, File dynRestartFile) {
169     this.potential = potentialEnergy;
170     this.orthogonalSpaceTempering = orthogonalSpaceTempering;
171     mcLogLevel = verbose ? Level.INFO : Level.FINE;
172     mdVerbosity = verbose ? MDVerbosity.QUIET : MDVerbosity.SILENT;
173     stepsPerMove = cycleLength;
174     totalSteps = dynamics.getNumSteps();
175 
176     ThermostatEnum thermostat = dynamics.thermostat;
177     if (!thermostat.equals(ThermostatEnum.ADIABATIC)) {
178       logger.info(format(" MC-OST MD moves will use the Adiabatic thermostat (%s ignored).", thermostat));
179       dynamics.setThermostat(ThermostatEnum.ADIABATIC);
180     }
181 
182     IntegratorEnum integrator = dynamics.integrator;
183     if (!integrator.reversible || !integrator.deterministic) {
184       logger.info(format(" MC-OST MD moves require a reversible, deterministic integrator (Verlet replaced %s).", integrator));
185       dynamics.setIntegrator(IntegratorEnum.VERLET);
186     }
187 
188     // Create the MC MD and Lambda moves.
189     boolean useOST = properties.getBoolean("mdmove-full", false);
190     Potential mdPotential = useOST ? orthogonalSpaceTempering : potential;
191     mdMove = new MDMove(molecularAssembly, mdPotential, listener, dynamics, stepsPerMove, dynRestartFile);
192     if (properties.containsKey("randomseed")) {
193       int randomSeed = properties.getInt("randomseed", 0);
194       logger.info(format(" Setting random seed for lambdaMove to %d ", randomSeed));
195       lambdaMove = new LambdaMove(randomSeed, orthogonalSpaceTempering);
196       setRandomSeed(randomSeed);
197     } else {
198       lambdaMove = new LambdaMove(orthogonalSpaceTempering);
199     }
200 
201     // Configure discrete Lambda.
202     boolean discreteLambda = orthogonalSpaceTempering.getHistogram().hd.discreteLambda;
203     if (discreteLambda) {
204       lambdaMove.setContinuous(false);
205       // Set the move size to the Lambda bin width.
206       lambdaMove.setMoveSize(orthogonalSpaceTempering.getHistogram().hd.lambdaBinWidth);
207     }
208 
209     // The OST class will handle adding the time dependent bias.
210     orthogonalSpaceTempering.setPropagateLambda(false);
211     int frequency = properties.getInt("mc-ost-bias-frequency", 1);
212     if (frequency <= 1) {
213       biasDepositionFrequency = 1;
214     } else {
215       biasDepositionFrequency = frequency;
216       logger.info(format(" MC-OST will deposit a bias only once per %d MC cycles.", biasDepositionFrequency));
217     }
218   }
219 
220   /**
221    * Returns the current value of lambda
222    *
223    * @return lambda
224    */
225   public double getLambda() {
226     return lambda;
227   }
228 
229   /**
230    * Calls on the OST method set lambda to update lambda to the current value in this class
231    *
232    * @param lambda a double.
233    */
234   public void setLambda(double lambda) {
235     this.lambda = lambda;
236     orthogonalSpaceTempering.setLambda(lambda);
237   }
238 
239   public MolecularDynamics getMD() {
240     return mdMove.getMD();
241   }
242 
243   /**
244    * The goal is to sample lambda and coordinates (X) simultaneously to converge the ensemble average
245    * dU/dL for every state (lambda) along the thermodynamic path. Note that the order of 1 and 2 below
246    * can be swapped (i.e. run MD and then change lambda). Here the order is random for each trial.
247    * <p>
248    * 1.) Randomly change the value of Lambda.
249    * <p>
250    * 2.) At the proposed lambda, run a am MD trajectory to "move" coordinates and dU/dL.
251    * <p>
252    * 3.) Accept / Reject the Lambda + MD move using the total Hamiltonian (Kinetic energy + OST energy).
253    * <p>
254    * 4.) Add to the bias.
255    */
256   public void sampleOneStep() {
257     // Validate the starting value of lambda.
258     lambda = orthogonalSpaceTempering.getLambda();
259     lambda = lambdaMove.validateLambda(lambda);
260     orthogonalSpaceTempering.setLambda(lambda);
261 
262     int n = potential.getNumberOfVariables();
263     double[] gradient = new double[n];
264     double[] currentCoordinates = new double[n];
265     double[] proposedCoordinates = new double[n];
266     long numMoves = totalSteps / stepsPerMove;
267     int acceptMCOST = 0;
268 
269     // Initialize the current coordinates.
270     potential.getCoordinates(currentCoordinates);
271 
272     // The Histogram maintains the time-dependent bias.
273     Histogram histogram = orthogonalSpaceTempering.getHistogram();
274 
275     // Compute the Total OST Energy as the sum of the Force Field Energy and Bias Energy.
276     double currentOSTEnergy = orthogonalSpaceTempering.energyAndGradient(currentCoordinates, gradient);
277 
278     // Retrieve the computed dU/dL, Force Field Energy and Bias Energy.
279     double currentdUdL = orthogonalSpaceTempering.getForceFielddEdL();
280     double currentForceFieldEnergy = orthogonalSpaceTempering.getForceFieldEnergy();
281     double currentBiasEnergy = orthogonalSpaceTempering.getBiasEnergy();
282 
283     // Loop over the number of requested MC Moves.
284     for (int imove = 0; imove < numMoves; imove++) {
285       long totalMoveTime = -nanoTime();
286 
287       // Change lambda before or after the MD trajectory.
288       boolean lambdaBeforeMD = random.nextBoolean();
289 
290       if (logger.isLoggable(Level.FINE)) {
291         if (equilibration) {
292           logger.fine(format("\n Equilibration Round %d", imove + 1));
293         } else {
294           String moveType = lambdaBeforeMD ? "Single-step lambda plus MD move." : "Single-step MD plus lambda move.";
295           logger.fine(format("\n MC-OST Round %d: %s", imove + 1, moveType));
296         }
297         logger.fine(format(" Starting force field energy for move %16.8f", currentForceFieldEnergy));
298       }
299 
300       // Set the current and proposed lambda.
301       double currentLambda = orthogonalSpaceTempering.getLambda();
302       double proposedLambda;
303       if (equilibration) {
304         // The proposed lambda equals the current lambda because its frozen for equilibration.
305         proposedLambda = currentLambda;
306         singleStepMD();
307       } else {
308         if (lambdaBeforeMD) {
309           // Change lambda before MD.
310           proposedLambda = singleStepLambda();
311           singleStepMD();
312         } else {
313           // Change lambda after MD.
314           singleStepMD();
315           proposedLambda = singleStepLambda();
316         }
317       }
318 
319       // Get the starting and final kinetic energy for the MD move.
320       double currentKineticEnergy = mdMove.getInitialKinetic();
321       double proposedKineticEnergy = mdMove.getKineticEnergy();
322 
323       // Get the new coordinates.
324       potential.getCoordinates(proposedCoordinates);
325 
326       // Compute the Total OST energy and gradient as the sum of the force field energy and bias energy.
327       long proposedOSTEnergyTime = -nanoTime();
328       double proposedOSTEnergy;
329       try {
330         proposedOSTEnergy = orthogonalSpaceTempering.energyAndGradient(proposedCoordinates, gradient);
331       } catch (EnergyException e) {
332         mdMove.revertMove();
333         if (!equilibration) {
334           lambdaMove.revertMove();
335           lambda = currentLambda;
336         }
337         logger.log(Level.INFO, " Unstable MD Move skipped.");
338         continue;
339       }
340       proposedOSTEnergyTime += nanoTime();
341 
342       if (logger.isLoggable(Level.FINE)) {
343         logger.fine(format(" Time to complete MD OST energy method call %6.3f", proposedOSTEnergyTime * NS2SEC));
344       }
345 
346       // Retrieve the proposed dU/dL, force field energy and bias energy.
347       double proposeddUdL = orthogonalSpaceTempering.getForceFielddEdL();
348       double proposedForceFieldEnergy = orthogonalSpaceTempering.getForceFieldEnergy();
349       double proposedBiasEnergy = orthogonalSpaceTempering.getBiasEnergy();
350 
351       // The Metropolis criteria is based on the sum of the OST Energy and Kinetic Energy.
352       double currentTotalEnergy = currentOSTEnergy + currentKineticEnergy;
353       double proposedTotalEnergy = proposedOSTEnergy + proposedKineticEnergy;
354 
355       // Log energy differences for the proposed move.
356       if (logger.isLoggable(mcLogLevel)) {
357         logger.log(mcLogLevel, format("\n  %8s %12s %12s %12s %12s", "", "Kinetic", "Potential", "Bias", "Total"));
358         logger.log(mcLogLevel, format("  Current  %12.4f %12.4f %12.4f %12.4f", currentKineticEnergy,
359                 currentForceFieldEnergy, currentBiasEnergy, currentTotalEnergy));
360         logger.log(mcLogLevel, format("  Proposed %12.4f %12.4f %12.4f %12.4f", proposedKineticEnergy,
361                 proposedForceFieldEnergy, proposedBiasEnergy, proposedTotalEnergy));
362         logger.log(mcLogLevel, format("  Delta    %12.4f %12.4f %12.4f %12.4f",
363             proposedKineticEnergy - currentKineticEnergy,
364             proposedForceFieldEnergy - currentForceFieldEnergy, proposedBiasEnergy - currentBiasEnergy,
365             proposedTotalEnergy - currentTotalEnergy));
366       }
367 
368       // Monitor energy conservation.
369       double energyChange = mdMove.getEnergyChange();
370       if (abs(energyChange) > ENERGY_CONSERVATION_TOLERANCE) {
371         mdMove.revertMove();
372         if (!equilibration) {
373           lambdaMove.revertMove();
374           lambda = currentLambda;
375         }
376         logger.warning(" MC Move skipped due to lack of MD energy conservation.");
377         continue;
378       }
379 
380       // Initialize the result string to reflect rejection.
381       String result = "    ";
382       // Energy change.
383       double dE = proposedTotalEnergy - currentTotalEnergy;
384       if (equilibration) {
385         if (orthogonalSpaceTempering.insideHardWallConstraint(proposedLambda, proposeddUdL)
386             && evaluateMove(currentTotalEnergy, proposedTotalEnergy)) {
387           // Accept Equilibration MD move.
388           result = " -> ";
389           acceptMCOST++;
390           currentOSTEnergy = proposedOSTEnergy;
391           currentdUdL = proposeddUdL;
392           currentForceFieldEnergy = proposedForceFieldEnergy;
393           currentBiasEnergy = proposedBiasEnergy;
394           arraycopy(proposedCoordinates, 0, currentCoordinates, 0, n);
395         } else {
396           mdMove.revertMove();
397         }
398       } else {
399         if (orthogonalSpaceTempering.insideHardWallConstraint(proposedLambda, proposeddUdL)
400             && evaluateMove(currentTotalEnergy, proposedTotalEnergy)) {
401           // Accept MD = Lambda move.
402           result = " -> ";
403           acceptMCOST++;
404           lambda = proposedLambda;
405           currentdUdL = proposeddUdL;
406           currentForceFieldEnergy = proposedForceFieldEnergy;
407           arraycopy(proposedCoordinates, 0, currentCoordinates, 0, n);
408         } else {
409           lambdaMove.revertMove();
410           lambda = currentLambda;
411           mdMove.revertMove();
412         }
413         if (imove % biasDepositionFrequency == 0) {
414           histogram.addBias(currentdUdL);
415           if (logger.isLoggable(Level.FINE)) {
416             logger.fine(format(" Added Bias at move %d: [ L=%5.3f, FL=%9.3f] ", imove + 1, lambda, currentdUdL));
417           }
418         } else {
419           logger.info(format(" Cycle %d: skipping bias deposition.", imove + 1));
420         }
421 
422         // Compute the updated OST bias.
423         currentBiasEnergy = histogram.computeBiasEnergy(lambda, currentdUdL);
424 
425         // Update the current OST Energy to be the sum of the current force field energy and updated OST Bias.
426         currentOSTEnergy = currentForceFieldEnergy + currentBiasEnergy;
427 
428         boolean snapShot = lambda >= orthogonalSpaceTempering.getLambdaWriteOut();
429         long mdMoveNum = (imove + 1) * stepsPerMove;
430         orthogonalSpaceTempering.getHistogram().ld.stepsTaken += stepsPerMove;
431         EnumSet<MDWriteAction> written = mdMove.writeFilesForStep(mdMoveNum, snapShot, true);
432         if (written.contains(MDWriteAction.RESTART)) {
433           orthogonalSpaceTempering.writeAdditionalRestartInfo(false);
434         }
435       }
436 
437       // Log the move.
438       totalMoveTime += nanoTime();
439       double percent = (acceptMCOST * 100.0) / (imove + 1);
440       logger.info(format(" %4d [ L=%5.3f, dU/dL=%8.3f]%s[ L=%5.3f, dU/dL=%8.3f] ΔE=%12.4f (%5.1f%%) in %6.3f sec.",
441           imove + 1, currentLambda, currentdUdL, result,
442           proposedLambda, proposeddUdL, dE, percent, totalMoveTime * NS2SEC));
443     }
444   }
445 
446   /**
447    * The goal is to sample lambda and coordinates (X) separately to converge the ensemble average
448    * dU/dL for every state (lambda) along the thermodynamic path.
449    * <p>
450    * 1.) At a fixed lambda, run a defined length MD trajectory to "move" coordinates and dU/dL.
451    * <p>
452    * 2.) Accept / Reject the MD move using the total Hamiltonian (Kinetic energy + OST energy).
453    * <p>
454    * 3.) Randomly change the value of Lambda.
455    * <p>
456    * 4.) Accept / Reject the Lambda move using the OST energy.
457    * <p>
458    * 5.) Add a hill to the histogram.
459    */
460   public void sampleTwoStep() {
461     // Validate the starting value of lambda.
462     lambda = orthogonalSpaceTempering.getLambda();
463     lambda = lambdaMove.validateLambda(lambda);
464     orthogonalSpaceTempering.setLambda(lambda);
465 
466     int n = potential.getNumberOfVariables();
467     double[] gradient = new double[n];
468     double[] currentCoordinates = new double[n];
469     double[] proposedCoordinates = new double[n];
470     long numMoves = totalSteps / stepsPerMove;
471     int acceptLambda = 0;
472     int acceptMD = 0;
473 
474     // Initialize the current coordinates.
475     potential.getCoordinates(currentCoordinates);
476 
477     // Update time dependent bias.
478     Histogram histogram = orthogonalSpaceTempering.getHistogram();
479 
480     // Compute the current OST potential energy.
481     double currentOSTEnergy = orthogonalSpaceTempering.energyAndGradient(currentCoordinates, gradient);
482 
483     // Collect the current dU/dL, Force Field Energy and Bias Energy.
484     double currentdUdL = orthogonalSpaceTempering.getForceFielddEdL();
485     double currentForceFieldEnergy = orthogonalSpaceTempering.getForceFieldEnergy();
486     double currentBiasEnergy = orthogonalSpaceTempering.getBiasEnergy();
487 
488     // Initialize MC move instances.
489     for (int imove = 0; imove < numMoves; imove++) {
490       long totalMoveTime = -nanoTime();
491       long mdMoveAndEvalTime = -nanoTime();
492 
493       if (logger.isLoggable(Level.FINE)) {
494         if (equilibration) {
495           logger.fine(format("\n MD Equilibration Round %d", imove + 1));
496         } else {
497           logger.fine(format("\n MC Orthogonal Space Sampling Round %d: Independent Steps", imove + 1));
498         }
499       }
500 
501       // Run MD in an approximate potential U* (U star) that does not include the OST bias.
502       long mdMoveTime = -nanoTime();
503       mdMove.move(mdVerbosity);
504       mdMoveTime += nanoTime();
505       logger.log(mcLogLevel, format("  Total time for MD move: %6.3f", mdMoveTime * NS2SEC));
506 
507       // Get the starting and final kinetic energy for the MD move.
508       double currentKineticEnergy = mdMove.getInitialKinetic();
509       double proposedKineticEnergy = mdMove.getKineticEnergy();
510 
511       // Get the new coordinates.
512       potential.getCoordinates(proposedCoordinates);
513 
514       // Compute the Total OST Energy as the sum of the Force Field Energy and Bias Energy.
515       long proposedOSTEnergyTime = -nanoTime();
516 
517       double proposedOSTEnergy;
518       try {
519         proposedOSTEnergy = orthogonalSpaceTempering.energyAndGradient(proposedCoordinates, gradient);
520       } catch (EnergyException e) {
521         mdMove.revertMove();
522         logger.log(Level.INFO, " Unstable MD Move skipped.");
523         continue;
524       }
525       proposedOSTEnergyTime += nanoTime();
526 
527       logger.fine(format("  Time to complete MD OST energy method call %6.3f", proposedOSTEnergyTime * NS2SEC));
528 
529       // Retrieve the proposed dU/dL, Force Field Energy and Bias Energy.
530       double proposeddUdL = orthogonalSpaceTempering.getForceFielddEdL();
531       double proposedForceFieldEnergy = orthogonalSpaceTempering.getForceFieldEnergy();
532       double proposedBiasEnergy = orthogonalSpaceTempering.getBiasEnergy();
533 
534       // The Metropolis criteria is based on the sum of the OST Energy and Kinetic Energy.
535       double currentTotalEnergy = currentOSTEnergy + currentKineticEnergy;
536       double proposedTotalEnergy = proposedOSTEnergy + proposedKineticEnergy;
537 
538       if (logger.isLoggable(mcLogLevel)) {
539         logger.log(mcLogLevel, format("\n  %8s %12s %12s %12s %12s", "", "Kinetic", "Potential", "Bias", "Total"));
540         logger.log(mcLogLevel, format("  Current  %12.4f %12.4f %12.4f %12.4f", currentKineticEnergy,
541             currentForceFieldEnergy, currentBiasEnergy, currentTotalEnergy));
542         logger.log(mcLogLevel, format("  Proposed %12.4f %12.4f %12.4f %12.4f", proposedKineticEnergy,
543             proposedForceFieldEnergy, proposedBiasEnergy, proposedTotalEnergy));
544         logger.log(mcLogLevel, format("  Delta    %12.4f %12.4f %12.4f %12.4f",
545             proposedKineticEnergy - currentKineticEnergy,
546             proposedForceFieldEnergy - currentForceFieldEnergy, proposedBiasEnergy - currentBiasEnergy,
547             proposedTotalEnergy - currentTotalEnergy));
548       }
549 
550       double energyChange = mdMove.getEnergyChange();
551       if (abs(energyChange) > ENERGY_CONSERVATION_TOLERANCE) {
552         mdMove.revertMove();
553         logger.warning(" MC Move skipped due to lack of MD energy conservation");
554         continue;
555       }
556 
557       String result = "    ";
558       double dE = proposedTotalEnergy - currentTotalEnergy;
559       if (orthogonalSpaceTempering.insideHardWallConstraint(orthogonalSpaceTempering.getLambda(),
560           proposeddUdL) && evaluateMove(currentTotalEnergy, proposedTotalEnergy)) {
561         // Accept MD move.
562         acceptMD++;
563         result = " -> ";
564         currentOSTEnergy = proposedOSTEnergy;
565         currentdUdL = proposeddUdL;
566         currentForceFieldEnergy = proposedForceFieldEnergy;
567         currentBiasEnergy = proposedBiasEnergy;
568         arraycopy(proposedCoordinates, 0, currentCoordinates, 0, n);
569       } else {
570         mdMove.revertMove();
571       }
572 
573       // Log the move.
574       mdMoveAndEvalTime += nanoTime();
575       double percent = (acceptMD * 100.0) / (imove + 1);
576       logger.info(format(" %4d MD [ L=%5.3f, dU/dL=%8.3f]%s[ L=%5.3f, dU/dL=%8.3f] ΔE=%12.4f (%5.1f%%) in %6.3f sec.", imove + 1,
577           lambda, currentdUdL, result, lambda, proposeddUdL, dE, percent, mdMoveAndEvalTime * NS2SEC));
578 
579       // During equilibration, do not change Lambda or contribute to the OST bias.
580       if (!equilibration) {
581         // Update Lambda.
582         long lambdaMoveTime = -nanoTime();
583         double currentLambda = orthogonalSpaceTempering.getLambda();
584         lambdaMove.move();
585         double proposedLambda = orthogonalSpaceTempering.getLambda();
586 
587         // Compute the Total OST Energy as the sum of the Force Field Energy and Bias Energy.
588         long proposedOSTEnergyTime2 = -nanoTime();
589         proposedOSTEnergy = orthogonalSpaceTempering.energyAndGradient(currentCoordinates, gradient);
590         proposedOSTEnergyTime2 += nanoTime();
591 
592         // Retrieve the proposed dU/dL, Force Field Energy and Bias Energy.
593         proposedForceFieldEnergy = orthogonalSpaceTempering.getForceFieldEnergy();
594         proposeddUdL = orthogonalSpaceTempering.getForceFielddEdL();
595 
596         if (logger.isLoggable(mcLogLevel)) {
597           logger.log(mcLogLevel, format("\n  Time to complete Lambda OST energy method call %6.3f ", proposedOSTEnergyTime2 * NS2SEC));
598           logger.log(mcLogLevel, format("  Current  OST     %12.3f at L=%5.3f.", currentOSTEnergy, currentLambda));
599           logger.log(mcLogLevel, format("  Proposed OST     %12.3f at L=%5.3f.", proposedOSTEnergy, proposedLambda));
600           logger.log(mcLogLevel, format("  MC Energy change: %12.3f (kcal/mol).", proposedOSTEnergy - currentOSTEnergy));
601         }
602 
603         dE = proposedOSTEnergy - currentOSTEnergy;
604         if (orthogonalSpaceTempering.insideHardWallConstraint(proposedLambda, proposeddUdL)
605             && evaluateMove(currentOSTEnergy, proposedOSTEnergy)) {
606           acceptLambda++;
607           result = " -> ";
608           currentForceFieldEnergy = proposedForceFieldEnergy;
609           currentdUdL = proposeddUdL;
610           lambda = proposedLambda;
611         } else {
612           result = "    ";
613           lambdaMove.revertMove();
614           lambda = currentLambda;
615         }
616 
617         lambdaMoveTime += nanoTime();
618         percent = (acceptLambda * 100.0) / (imove + 1);
619         logger.info(format("      L  [ L=%5.3f, dU/dL=%8.3f]%s[ L=%5.3f, dU/dL=%8.3f] ΔE=%12.4f (%5.1f%%) in %6.3f sec.",
620             currentLambda, currentdUdL, result, proposedLambda, proposeddUdL, dE, percent, lambdaMoveTime * NS2SEC));
621 
622         if (imove % biasDepositionFrequency == 0) {
623           histogram.addBias(currentdUdL);
624           logger.fine(format(" Added Bias at move %d: [ L=%5.3f, FL=%9.3f] ", imove + 1, lambda, currentdUdL));
625         } else {
626           logger.info(format(" Cycle %d: skipping bias deposition.", imove + 1));
627         }
628 
629         // Compute the updated OST bias.
630         currentBiasEnergy = histogram.computeBiasEnergy(lambda, currentdUdL);
631 
632         // Update the current OST Energy to be the sum of the current Force Field Energy and updated
633         // OST Bias.
634         currentOSTEnergy = currentForceFieldEnergy + currentBiasEnergy;
635 
636         if (automaticWriteouts) {
637           long mdMoveNum = (imove + 1) * stepsPerMove;
638           boolean trySnapshot = lambda >= orthogonalSpaceTempering.getLambdaWriteOut();
639           EnumSet<MDWriteAction> written = mdMove.writeFilesForStep(mdMoveNum, trySnapshot, true);
640           if (written.contains(MDWriteAction.RESTART)) {
641             orthogonalSpaceTempering.writeAdditionalRestartInfo(false);
642           }
643         }
644       }
645 
646       if (logger.isLoggable(Level.FINE)) {
647         totalMoveTime += nanoTime();
648         logger.fine(format(" Round complete in %6.3f sec.", totalMoveTime * NS2SEC));
649       }
650     }
651   }
652 
653   public void setAutomaticWriteouts(boolean autoWrite) {
654     automaticWriteouts = autoWrite;
655   }
656 
657   /**
658    * Sets the value of the boolean equilibration variables to true or false to either allow an
659    * equilibration step or skip it.
660    *
661    * @param equilibration a boolean.
662    */
663   public void setEquilibration(boolean equilibration) {
664     this.equilibration = equilibration;
665   }
666 
667   /**
668    * Calls on LambdaMove class method setLambdaStdDev to update the lambda standard deviation to the
669    * current value in this class
670    *
671    * @param stdDev a double.
672    */
673   public void setLambdaStdDev(double stdDev) {
674     if (!lambdaMove.isContinuous()) {
675       if (stdDev != orthogonalSpaceTempering.getHistogram().hd.lambdaBinWidth) {
676         logger.info(
677             format(" Requested Lambda step size change %6.4f is not equal to OST bin width %6.4f.",
678                 stdDev, orthogonalSpaceTempering.getHistogram().hd.lambdaBinWidth));
679         return;
680       }
681     }
682     lambdaMove.setMoveSize(stdDev);
683   }
684 
685   /**
686    * Sets the next number of steps MC-OST should run.
687    *
688    * @param totalSteps The length of the next MC-OST run.
689    */
690   public void setTotalSteps(long totalSteps) {
691     this.totalSteps = totalSteps;
692   }
693 
694   /**
695    * Propose a lambda move.
696    *
697    * @return The proposed lambda.
698    */
699   private double singleStepLambda() {
700     lambdaMove.move();
701     double proposedLambda = orthogonalSpaceTempering.getLambda();
702     logger.log(mcLogLevel, format(" Proposed lambda: %5.3f.", proposedLambda));
703     return proposedLambda;
704   }
705 
706   /**
707    * Run MD in an approximate potential U* (U star) that does not include the OST bias.
708    */
709   private void singleStepMD() {
710     long mdMoveTime = -nanoTime();
711     mdMove.move(mdVerbosity);
712     mdMoveTime += nanoTime();
713     logger.log(mcLogLevel, format(" Total time for MD move: %6.3f", mdMoveTime * NS2SEC));
714   }
715 
716   /**
717    * {@inheritDoc}
718    */
719   @Override
720   public void revertStep() {
721     throw new UnsupportedOperationException("Not supported yet.");
722   }
723 
724   /**
725    * {@inheritDoc}
726    */
727   @Override
728   protected double currentEnergy() {
729     throw new UnsupportedOperationException("Not supported yet.");
730   }
731 
732   /**
733    * {@inheritDoc}
734    */
735   @Override
736   protected void storeState() {
737     throw new UnsupportedOperationException("Not supported yet.");
738   }
739 }