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.dynamics;
39  
40  import static ffx.utilities.Constants.FSEC_TO_PSEC;
41  import static ffx.utilities.Constants.KCAL_TO_GRAM_ANG2_PER_PS2;
42  import static ffx.utilities.Constants.NS2SEC;
43  import static java.lang.String.format;
44  import static java.util.Arrays.fill;
45  
46  import edu.rit.pj.Comm;
47  import ffx.algorithms.AlgorithmListener;
48  import ffx.algorithms.Terminatable;
49  import ffx.algorithms.dynamics.integrators.BetterBeeman;
50  import ffx.algorithms.dynamics.integrators.Integrator;
51  import ffx.algorithms.dynamics.integrators.IntegratorEnum;
52  import ffx.algorithms.dynamics.integrators.Respa;
53  import ffx.algorithms.dynamics.integrators.Stochastic;
54  import ffx.algorithms.dynamics.integrators.VelocityVerlet;
55  import ffx.algorithms.dynamics.thermostats.Adiabatic;
56  import ffx.algorithms.dynamics.thermostats.Berendsen;
57  import ffx.algorithms.dynamics.thermostats.Bussi;
58  import ffx.algorithms.dynamics.thermostats.Thermostat;
59  import ffx.algorithms.dynamics.thermostats.ThermostatEnum;
60  import ffx.algorithms.thermodynamics.OrthogonalSpaceTempering;
61  import ffx.crystal.Crystal;
62  import ffx.numerics.Constraint;
63  import ffx.numerics.Potential;
64  import ffx.numerics.math.RunningStatistics;
65  import ffx.potential.openmm.OpenMMEnergy;
66  import ffx.potential.MolecularAssembly;
67  import ffx.potential.SystemState;
68  import ffx.potential.UnmodifiableState;
69  import ffx.potential.bonded.Atom;
70  import ffx.potential.bonded.LambdaInterface;
71  import ffx.potential.extended.ExtendedSystem;
72  import ffx.potential.parameters.ForceField;
73  import ffx.potential.parsers.DYNFilter;
74  import ffx.potential.parsers.PDBFilter;
75  import ffx.potential.parsers.XPHFilter;
76  import ffx.potential.parsers.XYZFilter;
77  import ffx.utilities.FileUtils;
78  import java.io.File;
79  import java.util.ArrayList;
80  import java.util.Arrays;
81  import java.util.EnumSet;
82  import java.util.List;
83  import java.util.logging.Level;
84  import java.util.logging.Logger;
85  import org.apache.commons.configuration2.CompositeConfiguration;
86  import org.apache.commons.io.FilenameUtils;
87  
88  /**
89   * Run NVE, NVT, or NPT molecular dynamics.
90   *
91   * @author Michael J. Schnieders
92   * @since 1.0
93   */
94  public class MolecularDynamics implements Runnable, Terminatable {
95  
96    private static final Logger logger = Logger.getLogger(MolecularDynamics.class.getName());
97  
98    /**
99     * The default interval (in picoseconds) between logging the current state.
100    */
101   private static final double DEFAULT_LOG_INTERVAL = 0.25;
102   /**
103    * The default interval (in picoseconds) between writing the current state to a DYN file.
104    */
105   private static final double DEFAULT_RESTART_INTERVAL = 1.0;
106   /**
107    * The default interval (in picoseconds) between writing the current coordinates.
108    */
109   private static final double DEFAULT_TRAJECTORY_INTERVAL = 10.0;
110   /** By default, wait 25000 nsec (25 usec) in between polling the dynamics thread. */
111   private static final int DEFAULT_DYNAMICS_SLEEP_TIME = 25000;
112 
113   /**
114    * MolecularAssembly to run dynamics on.
115    * <p>
116    * For dual or quad topology simulations, this is an array of length 2 or 4.
117    */
118   protected MolecularAssembly[] molecularAssembly;
119   /** Propagate dynamics on this potential surface. */
120   private final Potential potential;
121   /** An Algorithm Listener to send updates to the GUI. */
122   protected final AlgorithmListener algorithmListener;
123   /** Integrator instance. */
124   private final Integrator integrator;
125   /** Thermostat instance. */
126   private Thermostat thermostat;
127   /** Flag to indicate use of constant pressure. */
128   boolean constantPressure = false;
129   /**
130    * The Barostat in use, or null if constantPressure is false.
131    */
132   Barostat barostat = null;
133   /**
134    * Stores the current molecular dynamics state.
135    */
136   protected final SystemState state;
137   /**
138    * State of the system as of the last init call (and start of the last dynamics run).
139    */
140   protected UnmodifiableState initialState;
141   /**
142    * A copy of the current MD State. This can be used to revert to a previous state after a rejected
143    * MC Move.
144    */
145   private UnmodifiableState storedState;
146 
147   /**
148    * Temperature Stats for logging.
149    */
150   private final RunningStatistics temperatureStats = new RunningStatistics();
151   /**
152    * Potential Energy Stats for logging.
153    */
154   private final RunningStatistics potentialEnergyStats = new RunningStatistics();
155   /**
156    * Kinetic Energy Stats for logging.
157    */
158   private final RunningStatistics kineticEnergyStats = new RunningStatistics();
159   /**
160    * Total Energy Stats for logging.
161    */
162   private final RunningStatistics totalEnergyStats = new RunningStatistics();
163 
164   /**
165    * The time step (picoseconds).
166    * <p>
167    * The method <code>setTimeStep</code> should always be used to set this value.
168    */
169   protected double dt = 0.001;
170   /**
171    * The number of MD steps to take.
172    */
173   private long nSteps = 1000;
174   /**
175    * Target temperature. ToDo: use the Thermostat instance.
176    */
177   double targetTemperature = 298.15;
178   /**
179    * Flag to indicate velocities should be initialized.
180    */
181   protected boolean initVelocities = true;
182   /**
183    * Whether MD handles writing restart/trajectory files itself (true), or will be commanded by
184    * another class (false) to do it. The latter is true for MC-OST, for example.
185    */
186   protected boolean automaticWriteouts = true;
187   /**
188    * The total simulation time.
189    */
190   protected double totalSimTime = 0.0;
191 
192   /** Time between writing out restart/checkpoint files in picoseconds. */
193   protected double restartInterval = DEFAULT_RESTART_INTERVAL;
194   /** Time steps between writing out restart/checkpoint files. Set by the init method. */
195   protected int restartFrequency;
196   /** Time between appending to the trajectory file in picoseconds. */
197   private double trajectoryInterval = DEFAULT_TRAJECTORY_INTERVAL;
198   /** Time steps between adding a frame to the trajectory file. Set by the init method. */
199   protected int trajectoryFrequency;
200   /** Time between logging information to the screen in picoseconds. */
201   private double logInterval = DEFAULT_LOG_INTERVAL;
202   /**
203    * Time steps between logging information to the screen. Set by the init method.
204    */
205   protected int logFrequency;
206 
207   /** File type to use when saving files. */
208   protected String fileType = "XYZ";
209   /** Save snapshots in PDB format. */
210   boolean saveSnapshotAsPDB = true;
211   /** Dynamics restart file. */
212   File restartFile = null;
213   /** Flag to indicate loading of restart file. */
214   boolean loadRestart = false;
215   /** Filter to parse the dynamics restart file. */
216   DYNFilter dynFilter = null;
217   /** If asked to perform dynamics with a null dynamics file, write here. */
218   private File fallbackDynFile;
219 
220   /**
221    * Log basic information at this level.
222    */
223   Level basicLogging = Level.INFO;
224   /** Indicates how verbose MD should be. */
225   private MDVerbosity verbosityLevel = MDVerbosity.VERBOSE;
226 
227   /** Flag to indicate dynamics has been initialized. */
228   boolean initialized = false;
229   /** Flag to indicate MD should be terminated. */
230   private boolean terminate = false;
231   /** Flag to indicate a run has finished. */
232   protected boolean done;
233 
234   /** ESV System. */
235   private ExtendedSystem esvSystem;
236   /** ESV Stochastic integrator. */
237   private Stochastic esvIntegrator;
238   /** ESV Adiabatic thermostat. */
239   private Adiabatic esvThermostat;
240   /** Frequency to print ESV info. */
241   private int printEsvFrequency = -1;
242 
243   /**
244    * Constructor for MolecularDynamics.
245    *
246    * @param assembly a {@link ffx.potential.MolecularAssembly} object.
247    * @param potentialEnergy a {@link ffx.numerics.Potential} object.
248    * @param listener a {@link ffx.algorithms.AlgorithmListener} object.
249    * @param requestedThermostat a {@link ThermostatEnum} object.
250    * @param requestedIntegrator a {@link ffx.algorithms.dynamics.integrators.IntegratorEnum}
251    *     object.
252    */
253   public MolecularDynamics(MolecularAssembly assembly, Potential potentialEnergy,
254       AlgorithmListener listener, ThermostatEnum requestedThermostat,
255       IntegratorEnum requestedIntegrator) {
256     this.molecularAssembly = new MolecularAssembly[1];
257     this.molecularAssembly[0] = assembly;
258     this.algorithmListener = listener;
259     this.potential = potentialEnergy;
260     if (potentialEnergy instanceof Barostat) {
261       constantPressure = true;
262       barostat = (Barostat) potentialEnergy;
263     }
264 
265     int numberOfVariables = potentialEnergy.getNumberOfVariables();
266     state = new SystemState(numberOfVariables);
267     state.setMass(potential.getMass());
268 
269     CompositeConfiguration properties = molecularAssembly[0].getProperties();
270     // If an Integrator wasn't passed to the MD constructor, check for one specified as a property.
271     if (requestedIntegrator == null) {
272       String integrate = properties.getString("integrate", "VERLET").trim();
273       try {
274         integrate = integrate.toUpperCase().replace("-", "_");
275         requestedIntegrator = IntegratorEnum.valueOf(integrate.toUpperCase());
276       } catch (Exception e) {
277         requestedIntegrator = IntegratorEnum.VERLET;
278       }
279     }
280 
281     boolean oMMLogging = potential instanceof OpenMMEnergy;
282     List<Constraint> constraints = potentialEnergy.getConstraints();
283 
284     // Set the integrator.
285     integrator = switch (requestedIntegrator) {
286       case RESPA, MTS -> {
287         Respa respa = new Respa(state);
288         int in = molecularAssembly[0].getProperties().getInt("respa-dt", 4);
289         if (in < 2) {
290           in = 2;
291         }
292         if (!oMMLogging) {
293           respa.setInnerTimeSteps(in);
294         }
295         logger.log(Level.FINE, format(" Created a RESPA integrator with %d inner time steps.", in));
296         yield respa;
297       }
298       case STOCHASTIC, LANGEVIN -> {
299         double friction = properties.getDouble("friction", 91.0);
300         logger.log(Level.FINE, format(" Friction set at %.3f collisions/picosecond", friction));
301         Stochastic stochastic = new Stochastic(friction, state);
302         if (properties.containsKey("randomseed")) {
303           stochastic.setRandomSeed(properties.getInt("randomseed", 0));
304         }
305         // The stochastic dynamics integration procedure will thermostat
306         // the system. The ADIABTIC thermostat just serves to report the
307         // temperature and initialize velocities if necessary.
308         requestedThermostat = ThermostatEnum.ADIABATIC;
309         yield stochastic;
310       }
311       case BEEMAN -> new BetterBeeman(state);
312       default -> new VelocityVerlet(state);
313     };
314 
315     integrator.addConstraints(constraints);
316 
317     // If a Thermostat wasn't passed to the MD constructor, check for one specified as a property.
318     if (requestedThermostat == null) {
319       String thermo = properties.getString("thermostat", "Berendsen").trim();
320       try {
321         thermo = thermo.toUpperCase().replace("-", "_");
322         requestedThermostat = ThermostatEnum.valueOf(thermo.toUpperCase());
323       } catch (Exception e) {
324         requestedThermostat = ThermostatEnum.BERENDSEN;
325       }
326     }
327 
328     // Set the thermostat.
329     thermostat = switch (requestedThermostat) {
330       case BERENDSEN -> {
331         double tau = properties.getDouble("tau-temperature", 0.2);
332         yield new Berendsen(state, potentialEnergy.getVariableTypes(), targetTemperature, tau,
333             constraints);
334       }
335       case BUSSI -> {
336         double tau = properties.getDouble("tau-temperature", 0.2);
337         Bussi bussi = new Bussi(state, potentialEnergy.getVariableTypes(), targetTemperature, tau,
338             constraints);
339         if (properties.containsKey("randomseed")) {
340           bussi.setRandomSeed(properties.getInt("randomseed", 0));
341         }
342         yield bussi;
343       }
344       default -> new Adiabatic(state, potentialEnergy.getVariableTypes(), constraints);
345     };
346 
347     if (properties.containsKey("randomseed")) {
348       thermostat.setRandomSeed(properties.getInt("randomseed", 0));
349     }
350 
351     // For Stochastic dynamics, center of mass motion will not be removed.
352     if (integrator instanceof Stochastic) {
353       thermostat.setRemoveCenterOfMassMotion(false);
354     }
355 
356     done = true;
357     fallbackDynFile = defaultFallbackDyn(assembly);
358   }
359 
360   /**
361    * Method that determines whether a dynamics is done by the java implementation native to ffx or
362    * the OpenMM implementation
363    *
364    * @param assembly a {@link ffx.potential.MolecularAssembly} object.
365    * @param potentialEnergy a {@link ffx.numerics.Potential} object.
366    * @param listener a {@link ffx.algorithms.AlgorithmListener} object.
367    * @param requestedThermostat a {@link ThermostatEnum} object.
368    * @param requestedIntegrator a {@link ffx.algorithms.dynamics.integrators.IntegratorEnum}
369    *     object.
370    * @return a {@link MolecularDynamics} object.
371    */
372   public static MolecularDynamics dynamicsFactory(MolecularAssembly assembly,
373       Potential potentialEnergy, AlgorithmListener listener, ThermostatEnum requestedThermostat,
374       IntegratorEnum requestedIntegrator) {
375 
376     return dynamicsFactory(assembly, potentialEnergy, listener, requestedThermostat,
377         requestedIntegrator, defaultEngine(assembly, potentialEnergy));
378   }
379 
380   /**
381    * dynamicsFactory.
382    *
383    * @param assembly a {@link ffx.potential.MolecularAssembly} object.
384    * @param potentialEnergy a {@link ffx.numerics.Potential} object.
385    * @param listener a {@link ffx.algorithms.AlgorithmListener} object.
386    * @param requestedThermostat a {@link ThermostatEnum} object.
387    * @param requestedIntegrator a {@link ffx.algorithms.dynamics.integrators.IntegratorEnum}
388    *     object.
389    * @param engine a {@link MDEngine} object.
390    * @return a {@link MolecularDynamics} object.
391    */
392   public static MolecularDynamics dynamicsFactory(MolecularAssembly assembly,
393       Potential potentialEnergy, AlgorithmListener listener, ThermostatEnum requestedThermostat,
394       IntegratorEnum requestedIntegrator, MDEngine engine) {
395 
396     switch (engine) {
397       case OPENMM, OMM -> {
398         // TODO: Replace this with calls to the leaves of a proper tree structure.
399         // Unfortunately, neither Java, nor Apache Commons, nor Guava has an arbitrary tree
400         // implementing Collection.
401         // Nor does javax.swing have a quick "get me the leaves" method that I was able to find.
402         boolean ommLeaves = potentialEnergy.getUnderlyingPotentials().stream()
403             .anyMatch((Potential p) -> p instanceof OpenMMEnergy);
404         ommLeaves = ommLeaves || potentialEnergy instanceof OpenMMEnergy;
405         if (ommLeaves) {
406           return new MolecularDynamicsOpenMM(assembly, potentialEnergy, listener,
407               requestedThermostat, requestedIntegrator);
408         } else {
409           throw new IllegalArgumentException(format(
410               " Requested OpenMM engine %s, but at least one leaf of the potential %s is not an OpenMM force field!",
411               engine, potentialEnergy));
412         }
413       }
414       default -> {
415         return new MolecularDynamics(assembly, potentialEnergy, listener, requestedThermostat,
416             requestedIntegrator);
417       }
418     }
419   }
420 
421   private static MDEngine defaultEngine(MolecularAssembly molecularAssembly,
422       Potential potentialEnergy) {
423     CompositeConfiguration properties = molecularAssembly.getProperties();
424     String mdEngine = properties.getString("MD-engine");
425     if (mdEngine != null) {
426       if (mdEngine.equalsIgnoreCase("OMM")) {
427         logger.info(" Creating OpenMM Dynamics Object");
428         return MDEngine.OPENMM;
429       } else {
430         logger.info(" Creating FFX Dynamics Object");
431         return MDEngine.FFX;
432       }
433     } else {
434       // TODO: Replace this with a better check.
435       boolean ommLeaves = potentialEnergy.getUnderlyingPotentials().stream()
436           .anyMatch((Potential p) -> p instanceof OpenMMEnergy);
437       ommLeaves = ommLeaves || potentialEnergy instanceof OpenMMEnergy;
438       if (ommLeaves) {
439         return MDEngine.OPENMM;
440       } else {
441         return MDEngine.FFX;
442       }
443     }
444   }
445 
446   /**
447    * Creates the "default" fallback dynamics file object (does not create actual file!)
448    *
449    * @param assembly First assembly.
450    * @return Default fallback file.
451    */
452   private static File defaultFallbackDyn(MolecularAssembly assembly) {
453     String firstFileName = FilenameUtils.removeExtension(assembly.getFile().getAbsolutePath());
454     return new File(firstFileName + ".dyn");
455   }
456 
457   /**
458    * Adds a MolecularAssembly to be tracked by this MolecularDynamics. Note: does not affect the
459    * underlying Potential.
460    *
461    * @param assembly A MolecularAssembly to be tracked.
462    */
463   public void addAssembly(MolecularAssembly assembly) {
464     molecularAssembly = Arrays.copyOf(molecularAssembly, molecularAssembly.length + 1);
465     molecularAssembly[molecularAssembly.length - 1] = assembly;
466   }
467 
468   /**
469    * attachExtendedSystem.
470    *
471    * @param system a {@link ffx.potential.extended.ExtendedSystem} object.
472    */
473   public void attachExtendedSystem(ExtendedSystem system, double reportFreq) {
474     if (esvSystem != null) {
475       logger.warning("An ExtendedSystem is already attached to this MD!");
476     }
477     esvSystem = system;
478     SystemState esvState = esvSystem.getState();
479     this.esvIntegrator = new Stochastic(esvSystem.getThetaFriction(), esvState);
480     if(!esvSystem.getConstraints().isEmpty()){
481       esvIntegrator.addConstraints(esvSystem.getConstraints());
482     }
483     this.esvThermostat = new Adiabatic(esvState, potential.getVariableTypes());
484     printEsvFrequency = intervalToFreq(reportFreq, "Reporting (logging) interval");
485     logger.info(
486         format("  Attached extended system (%s) to molecular dynamics.", esvSystem.toString()));
487     logger.info(format("  Extended System Theta Friction: %f", esvSystem.getThetaFriction()));
488     logger.info(format("  Extended System Theta Mass: %f", esvSystem.getThetaMass()));
489     logger.info(format("  Extended System Lambda Print Frequency: %d (fsec)", printEsvFrequency));
490   }
491 
492   /**
493    * Blocking molecular dynamics. When this method returns, the MD run is done.
494    *
495    * @param nSteps Number of MD steps
496    * @param timeStep Time step in femtoseconds
497    * @param loggingInterval Interval between printing/logging information in picoseconds.
498    * @param trajectoryInterval Interval between adding a frame to the trajectory file in
499    *     picoseconds.
500    * @param temperature Temperature in Kelvins.
501    * @param initVelocities Initialize new velocities from a Maxwell-Boltzmann distribution.
502    * @param fileType XYZ or ARC to save to .arc, PDB for .pdb files
503    * @param restartInterval Interval between writing new restart files in picoseconds.
504    * @param dyn A {@link java.io.File} object to write the restart file to.
505    */
506   public void dynamic(final long nSteps, final double timeStep, final double loggingInterval,
507       final double trajectoryInterval, final double temperature, final boolean initVelocities,
508       String fileType, double restartInterval, final File dyn) {
509     this.fileType = fileType;
510     setRestartFrequency(restartInterval);
511     dynamic(nSteps, timeStep, loggingInterval, trajectoryInterval, temperature, initVelocities, dyn);
512   }
513 
514   /**
515    * Blocking molecular dynamics. When this method returns, the MD run is done.
516    *
517    * @param nSteps Number of MD steps
518    * @param timeStep Time step (fsec)
519    * @param loggingInterval Interval between printing/logging information in picoseconds.
520    * @param trajectoryInterval Interval between adding a frame to the trajectory file in
521    *     picoseconds.
522    * @param temperature Temperature in Kelvins.
523    * @param initVelocities Initialize new velocities from a Maxwell-Boltzmann distribution.
524    * @param dyn A {@link java.io.File} object to write the restart file to.
525    */
526   public void dynamic(final long nSteps, final double timeStep, final double loggingInterval,
527       final double trajectoryInterval, final double temperature, final boolean initVelocities,
528       final File dyn) {
529     // Return if already running;
530     // Could happen if two threads call dynamic on the same MolecularDynamics instance.
531     if (!done) {
532       logger.warning(" Programming error - a thread invoked dynamic when it was already running.");
533       return;
534     }
535 
536     init(nSteps, timeStep, loggingInterval, trajectoryInterval, fileType, restartInterval,
537         temperature, initVelocities, dyn);
538 
539     Thread dynamicThread = new Thread(this);
540     dynamicThread.start();
541     synchronized (this) {
542       try {
543         while (dynamicThread.isAlive()) {
544           wait(0, DEFAULT_DYNAMICS_SLEEP_TIME);
545         }
546       } catch (InterruptedException e) {
547         String message = " Molecular dynamics interrupted.";
548         logger.log(Level.WARNING, message, e);
549       }
550     }
551     if (!verbosityLevel.isQuiet()) {
552       logger.info(" Done with an MD round.");
553     }
554   }
555 
556   /**
557    * Returns the MolecularAssembly array.
558    *
559    * @return A copy of the MolecularAssembly array.
560    */
561   public MolecularAssembly[] getAssemblyArray() {
562     return molecularAssembly.clone();
563   }
564 
565   /**
566    * Returns the associated dynamics file.
567    *
568    * @return Dynamics restart File.
569    */
570   public File getDynFile() {
571     return restartFile;
572   }
573 
574   /**
575    * Gets the kinetic energy at the start of the last dynamics run.
576    *
577    * @return Kinetic energy at the start of the run.
578    */
579   public double getInitialKineticEnergy() {
580     return initialState.kineticEnergy();
581   }
582 
583   /**
584    * Gets the potential energy at the start of the last dynamics run.
585    *
586    * @return potential energy at the start of the run.
587    */
588   public double getInitialPotentialEnergy() {
589     return initialState.potentialEnergy();
590   }
591 
592   /**
593    * Gets the temperature at the start of the last dynamics run.
594    *
595    * @return temperature at the start of the run.
596    */
597   public double getInitialTemperature() {
598     return initialState.temperature();
599   }
600 
601   /**
602    * Gets the total energy at the start of the last dynamics run.
603    *
604    * @return total energy at the start of the run.
605    */
606   public double getInitialTotalEnergy() {
607     return initialState.getTotalEnergy();
608   }
609 
610   /**
611    * getIntervalSteps.
612    *
613    * @return Always 1 for this implementation.
614    */
615   public int getIntervalSteps() {
616     return 1;
617   }
618 
619   /**
620    * No-op; FFX does not need to occasionally return information from FFX.
621    *
622    * @param intervalSteps Ignored.
623    */
624   public void setIntervalSteps(int intervalSteps) {
625     // Not meaningful for FFX MD.
626   }
627 
628   /**
629    * Get the system kinetic energy.
630    *
631    * @return kinetic energy.
632    */
633   public double getKineticEnergy() {
634     return state.getKineticEnergy();
635   }
636 
637   /**
638    * Get the system potential energy.
639    *
640    * @return potential energy.
641    */
642   public double getPotentialEnergy() {
643     return state.getPotentialEnergy();
644   }
645 
646   /**
647    * Get the current temperature of the system
648    *
649    * @return currentTemperature
650    */
651   public double getTemperature() {
652     return state.getTemperature();
653   }
654 
655   /**
656    * Getter for the field <code>thermostat</code>.
657    *
658    * @return a {@link ffx.algorithms.dynamics.thermostats.Thermostat} object.
659    */
660   public Thermostat getThermostat() {
661     return thermostat;
662   }
663 
664   /**
665    * Setter for the field <code>thermostat</code>.
666    *
667    * @param thermostat a {@link ffx.algorithms.dynamics.thermostats.Thermostat} object.
668    */
669   public void setThermostat(Thermostat thermostat) {
670     this.thermostat = thermostat;
671   }
672 
673   /**
674    * getTimeStep.
675    *
676    * @return Time step in picoseconds.
677    */
678   public double getTimeStep() {
679     return dt;
680   }
681 
682   /**
683    * Sets the time step and resets frequencies based on stored intervals.
684    *
685    * @param step Time step in femtoseconds.
686    */
687   private void setTimeStep(double step) {
688     dt = step * FSEC_TO_PSEC;
689     // Reset frequencies to be consistent with new time step.
690     setRestartFrequency(restartInterval);
691     setLoggingFrequency(logInterval);
692     setTrajectoryFrequency(trajectoryInterval);
693   }
694 
695   /**
696    * Get the total system energy (kinetic plus potential).
697    *
698    * @return total energy.
699    */
700   public double getTotalEnergy() {
701     return state.getTotalEnergy();
702   }
703 
704   public MDVerbosity getVerbosityLevel() {
705     return verbosityLevel;
706   }
707 
708   public void setVerbosityLevel(MDVerbosity level) {
709     verbosityLevel = level;
710     if (level == MDVerbosity.SILENT) {
711       basicLogging = Level.FINE;
712     } else {
713       basicLogging = Level.INFO;
714     }
715   }
716 
717   /**
718    * init
719    *
720    * @param nSteps Number of MD steps
721    * @param timeStep Time step in femtoseconds
722    * @param loggingInterval Interval between printing/logging information in picoseconds.
723    * @param trajectoryInterval Interval between adding a frame to the trajectory file in
724    *     picoseconds.
725    * @param fileType XYZ or ARC to save to .arc, PDB for .pdb files
726    * @param restartInterval Interval between writing new restart files in picoseconds.
727    * @param temperature Temperature in Kelvins.
728    * @param initVelocities Initialize new velocities from a Maxwell-Boltzmann distribution.
729    * @param dyn A {@link java.io.File} object to write the restart file to.
730    */
731   public void init(final long nSteps, final double timeStep, final double loggingInterval,
732       final double trajectoryInterval, final String fileType, final double restartInterval,
733       final double temperature, final boolean initVelocities, final File dyn) {
734 
735     // Return if already running.
736     if (!done) {
737       logger.warning(
738           " Programming error - attempt to modify parameters of a running MolecularDynamics instance.");
739       return;
740     }
741 
742     if (integrator instanceof Stochastic) {
743       if (constantPressure) {
744         logger.log(basicLogging, "\n Stochastic dynamics in the NPT ensemble");
745       } else {
746         logger.log(basicLogging, "\n Stochastic dynamics in the NVT ensemble");
747       }
748     } else if (!(thermostat instanceof Adiabatic)) {
749       if (constantPressure) {
750         logger.log(basicLogging, "\n Molecular dynamics in the NPT ensemble");
751       } else {
752         logger.log(basicLogging, "\n Molecular dynamics in the NVT ensemble");
753       }
754     } else {
755       if (constantPressure) {
756         logger.severe("\n NPT Molecular dynamics requires a thermostat");
757       } else {
758         logger.log(basicLogging, "\n Molecular dynamics in the NVE ensemble");
759       }
760     }
761 
762     this.nSteps = nSteps;
763     totalSimTime = 0.0;
764 
765     // Convert the time step from femtoseconds to picoseconds.
766     setTimeStep(timeStep);
767 
768     // Convert intervals in ps to frequencies in time steps.
769     setLoggingFrequency(loggingInterval);
770     setTrajectoryFrequency(trajectoryInterval);
771     setRestartFrequency(restartInterval);
772 
773     checkFrequency("Reporting           ", logFrequency);
774     if (automaticWriteouts) {
775       // Only sanity check these values if MD is doing this itself.
776       checkFrequency("Trajectory Write-Out", trajectoryFrequency);
777       checkFrequency("Restart             ", restartFrequency);
778     }
779 
780     // Set snapshot file type.
781     saveSnapshotAsPDB = true;
782     if (fileType.equalsIgnoreCase("XYZ") || fileType.equalsIgnoreCase("ARC")) {
783       saveSnapshotAsPDB = false;
784     } else if (!fileType.equalsIgnoreCase("PDB")) {
785       logger.warning("Snapshot file type unrecognized; saving snapshots as PDB.\n");
786     }
787 
788     setArchiveFile();
789 
790     this.restartFile = (dyn == null) ? fallbackDynFile : dyn;
791     loadRestart = restartFile.exists() && !initialized;
792 
793     if (dynFilter == null) {
794       dynFilter = new DYNFilter(molecularAssembly[0].getName());
795     }
796 
797     this.targetTemperature = temperature;
798     this.initVelocities = initVelocities;
799     done = false;
800 
801     if (loadRestart) {
802       logger.info("  Continuing from " + restartFile.getAbsolutePath());
803     }
804 
805     if (!verbosityLevel.isQuiet()) {
806       logger.info(format("  Integrator:          %15s", integrator.toString()));
807       logger.info(format("  Thermostat:          %15s", thermostat.toString()));
808       logger.info(format("  Number of steps:     %8d", nSteps));
809       logger.info(format("  Time step:           %8.3f (fsec)", timeStep));
810       logger.info(format("  Print interval:      %8.3f (psec)", loggingInterval));
811       logger.info(format("  Save interval:       %8.3f (psec)", trajectoryInterval));
812       if (molecularAssembly.length > 1) {
813         for (int i = 0; i < molecularAssembly.length; i++) {
814           File archiveFile = molecularAssembly[i].getArchiveFile();
815           logger.info(format("  Archive file %3d: %s", (i + 1), archiveFile.getAbsolutePath()));
816         }
817       } else {
818         logger.info(format("  Archive file:     %s",
819             molecularAssembly[0].getArchiveFile().getAbsolutePath()));
820       }
821       logger.info(format("  Restart file:     %s", restartFile.getAbsolutePath()));
822     }
823   }
824 
825   /**
826    * A version of init with the original method header. Redirects to the new method with default
827    * values for added parameters. Needed by (at least) ReplicaExchange, which calls this directly.
828    *
829    * @param nSteps Number of MD steps
830    * @param timeStep Time step in femtoseconds
831    * @param loggingInterval Interval between printing/logging information in picoseconds.
832    * @param trajectoryInterval Interval between adding a frame to the trajectory file in
833    *     picoseconds.
834    * @param temperature Temperature in Kelvins.
835    * @param initVelocities Initialize new velocities from a Maxwell-Boltzmann distribution.
836    * @param dyn A {@link java.io.File} object to write the restart file to.
837    */
838   public void init(final long nSteps, final double timeStep, final double loggingInterval,
839       final double trajectoryInterval, final double temperature, final boolean initVelocities,
840       final File dyn) {
841     init(nSteps, timeStep, loggingInterval, trajectoryInterval, "XYZ", restartInterval, temperature,
842         initVelocities, dyn);
843   }
844 
845   /** {@inheritDoc} */
846   @Override
847   public void run() {
848     try {
849       preRunOps();
850     } catch (IllegalStateException ise) {
851       return;
852     }
853     initializeEnergies();
854     postInitEnergies();
855     mainLoop();
856     postRun();
857   }
858 
859   /**
860    * Enable or disable automatic writeout of trajectory snapshots and restart files.
861    *
862    * <p>Primarily intended for use with MC-OST, which handles the timing itself.
863    *
864    * @param autoWriteout Whether to automatically write out trajectory and restart files.
865    */
866   public void setAutomaticWriteouts(boolean autoWriteout) {
867     this.automaticWriteouts = autoWriteout;
868   }
869 
870   /**
871    * Sets the "fallback" .dyn file to write to if none is passed to the dynamic method.
872    *
873    * @param fallback Fallback dynamics restart file.
874    */
875   public void setFallbackDynFile(File fallback) {
876     this.fallbackDynFile = fallback;
877   }
878 
879   /**
880    * Method to set file type from groovy scripts.
881    *
882    * @param fileType the type of snapshot files to write.
883    */
884   public void setFileType(String fileType) {
885     this.fileType = fileType;
886   }
887 
888   /**
889    * Not meaningful for FFX dynamics (no need to obtain velocities/accelerations from a different
890    * program, especially one running on a GPU). Is a no-op.
891    *
892    * @param obtainVA Not meaningful for this implementation.
893    */
894   public void setObtainVelAcc(boolean obtainVA) {
895     // Not meaningful for FFX dynamics.
896   }
897 
898   /**
899    * Method to set the Restart Frequency.
900    *
901    * @param restartInterval Time in ps between writing restart files.
902    * @throws java.lang.IllegalArgumentException If restart frequency is not a positive number
903    */
904   public void setRestartFrequency(double restartInterval) throws IllegalArgumentException {
905     restartFrequency = intervalToFreq(restartInterval, "Restart interval");
906     this.restartInterval = restartInterval;
907   }
908 
909   /**
910    * Set the archive file for each MolecularAssembly.
911    *
912    * @param archiveFiles An array of archive files.
913    */
914   public void setArchiveFiles(File[] archiveFiles) {
915     logger.info(" Setting archive files:\n " + Arrays.toString(archiveFiles));
916     int n = molecularAssembly.length;
917     assert archiveFiles.length == n;
918     for (int i = 0; i < n; i++) {
919       molecularAssembly[i].setArchiveFile(archiveFiles[i]);
920     }
921   }
922 
923   /** {@inheritDoc} */
924   @Override
925   public void terminate() {
926     terminate = true;
927     while (!done) {
928       synchronized (this) {
929         try {
930           wait(1);
931         } catch (Exception e) {
932           logger.log(Level.WARNING, " Exception terminating dynamics.\n", e);
933         }
934       }
935     }
936   }
937 
938   /**
939    * Write restart and trajectory files if the provided step matches the frequency and that file type
940    * is requested.
941    *
942    * @param step Step to write files (if any) for.
943    * @param trySnapshot If false, do not write snapshot even if the time step is correct.
944    * @param tryRestart If false, do not write a restart file even if the time step is correct.
945    * @return EnumSet of actions taken by this method.
946    */
947   public EnumSet<MDWriteAction> writeFilesForStep(long step, boolean trySnapshot,
948       boolean tryRestart) {
949     return writeFilesForStep(step, trySnapshot, tryRestart, null);
950   }
951 
952   /**
953    * Write restart and trajectory files if the provided step matches the frequency and that file type
954    * is requested.
955    *
956    * @param step Step to write files (if any) for.
957    * @param trySnapshot If false, do not write snapshot even if the time step is correct.
958    * @param tryRestart If false, do not write a restart file even if the time step is correct.
959    * @param extraLines Additional lines to append into the comments section of the snapshot (or
960    *     null).
961    * @return EnumSet of actions taken by this method.
962    */
963   public EnumSet<MDWriteAction> writeFilesForStep(long step, boolean trySnapshot, boolean tryRestart,
964       String[] extraLines) {
965     List<String> linesList =
966         (extraLines == null) ? new ArrayList<>() : new ArrayList<>(Arrays.asList(extraLines));
967 
968     if (potential instanceof LambdaInterface) {
969       String lamString = format("Lambda: %.8f", ((LambdaInterface) potential).getLambda());
970       linesList.add(lamString);
971     }
972 
973     String tempString = format("Temp: %.2f", thermostat.getTargetTemperature());
974     linesList.add(tempString);
975 
976     String timeString = format(" Time: %7.3e", totalSimTime);
977     linesList.add(timeString);
978 
979     if (esvSystem != null) {
980       String pHString = format("pH: %.2f", esvSystem.getConstantPh());
981       linesList.add(pHString);
982     }
983 
984     Comm world = Comm.world();
985     if (world != null && world.size() > 1) {
986       String rankString = format("Rank: %d", world.rank());
987       linesList.add(rankString);
988     }
989 
990     String[] allLines = new String[linesList.size()];
991     allLines = linesList.toArray(allLines);
992 
993     EnumSet<MDWriteAction> written = EnumSet.noneOf(MDWriteAction.class);
994     if (step != 0) {
995       // Write out snapshots in selected format every saveSnapshotFrequency steps.
996       if (trySnapshot && trajectoryFrequency > 0 && step % trajectoryFrequency == 0) {
997 
998         // Log stats.
999         logger.info(format("\n Average Values for the Last %d Out of %d Dynamics Steps\n",
1000             trajectoryFrequency, step));
1001         logger.info(format("  Simulation Time  %16.4f Picosecond", step * dt));
1002         logger.info(format("  Total Energy     %16.4f Kcal/mole   (+/-%9.4f)", totalEnergyStats.getMean(),
1003                 totalEnergyStats.getStandardDeviation()));
1004         logger.info(format("  Potential Energy %16.4f Kcal/mole   (+/-%9.4f)",
1005             potentialEnergyStats.getMean(), potentialEnergyStats.getStandardDeviation()));
1006         logger.info(format("  Kinetic Energy   %16.4f Kcal/mole   (+/-%9.4f)", kineticEnergyStats.getMean(),
1007                 kineticEnergyStats.getStandardDeviation()));
1008         logger.info(format("  Temperature      %16.4f Kelvin      (+/-%9.4f)\n", temperatureStats.getMean(),
1009                 temperatureStats.getStandardDeviation()));
1010         totalEnergyStats.reset();
1011         potentialEnergyStats.reset();
1012         kineticEnergyStats.reset();
1013         temperatureStats.reset();
1014 
1015         if (esvSystem != null) {
1016           for (Atom atom : molecularAssembly[0].getAtomList()) {
1017             int atomIndex = atom.getIndex() - 1;
1018             atom.setOccupancy(esvSystem.getTitrationLambda(atomIndex));
1019             atom.setTempFactor(esvSystem.getTautomerLambda(atomIndex));
1020           }
1021         }
1022         appendSnapshot(allLines);
1023         written.add(MDWriteAction.SNAPSHOT);
1024       }
1025 
1026       // Write out restart files every saveRestartFileFrequency steps.
1027       if (tryRestart && restartFrequency > 0 && step % restartFrequency == 0) {
1028         writeRestart();
1029         written.add(MDWriteAction.RESTART);
1030       }
1031     }
1032     return written;
1033   }
1034 
1035   /** Write out a restart file. */
1036   public void writeRestart() {
1037     String dynName = FileUtils.relativePathTo(restartFile).toString();
1038     double[] x = state.x();
1039     double[] v = state.v();
1040     double[] a = state.a();
1041     double[] aPrevious = state.aPrevious();
1042     if (dynFilter.writeDYN(restartFile, molecularAssembly[0].getCrystal(), x, v, a, aPrevious)) {
1043       logger.log(basicLogging, format(" Wrote dynamics restart to:  %s.", dynName));
1044     } else {
1045       logger.log(basicLogging, format(" Writing dynamics restart failed:  %s.", dynName));
1046     }
1047     if (esvSystem != null) {
1048       esvSystem.writeRestart();
1049       esvSystem.writeLambdaHistogram(false);
1050     }
1051     potential.writeAdditionalRestartInfo(true);
1052   }
1053 
1054   /**
1055    * Set the archive file for each MolecularAssembly (if not already set).
1056    */
1057   private void setArchiveFile() {
1058     for (MolecularAssembly assembly : molecularAssembly) {
1059       File file = assembly.getFile();
1060       String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
1061       File archiveFile = assembly.getArchiveFile();
1062       if (archiveFile == null) {
1063         archiveFile = new File(filename + ".arc");
1064         assembly.setArchiveFile(XYZFilter.version(archiveFile));
1065       }
1066     }
1067   }
1068 
1069   /**
1070    * Converts an interval in ps to a frequency in time steps.
1071    *
1072    * @param interval Interval between events in ps
1073    * @param describe Description of the event.
1074    * @return Frequency of event in time steps per event.
1075    * @throws IllegalArgumentException If interval is not a positive finite value.
1076    */
1077   private int intervalToFreq(double interval, String describe) throws IllegalArgumentException {
1078     if (!Double.isFinite(interval) || interval <= 0) {
1079       throw new IllegalArgumentException(
1080           format(" %s must be " + "positive finite value in ps, was %10.4g", describe, interval));
1081     }
1082     if (interval >= dt) {
1083       return (int) (interval / dt);
1084     } else {
1085       logger.warning(format(" Specified %s of %.6f ps < time step %.6f ps; "
1086           + "interval is set to once per time step!", describe, interval, dt));
1087       return 1;
1088     }
1089   }
1090 
1091   /**
1092    * Method to set the logging/printing frequency.
1093    *
1094    * @param logInterval Time in ps between logging information to the screen.
1095    */
1096   private void setLoggingFrequency(double logInterval) {
1097     logFrequency = intervalToFreq(logInterval, "Reporting (logging) interval");
1098     this.logInterval = logInterval;
1099   }
1100 
1101   /**
1102    * Method to set the frequency of appending snapshots to the trajectory file.
1103    *
1104    * @param snapshotInterval Time in ps between appending snapshots to the trajectory file.
1105    */
1106   private void setTrajectoryFrequency(double snapshotInterval) {
1107     trajectoryFrequency = intervalToFreq(snapshotInterval, "Trajectory writeout interval");
1108     this.trajectoryInterval = snapshotInterval;
1109   }
1110 
1111   /** Performs basic pre-MD operations such as loading the restart file. */
1112   void preRunOps() throws IllegalStateException {
1113     done = false;
1114     terminate = false;
1115 
1116     // Set the target temperature.
1117     thermostat.setTargetTemperature(targetTemperature);
1118     boolean quiet = verbosityLevel.isQuiet();
1119     thermostat.setQuiet(quiet);
1120     if (integrator instanceof Stochastic stochastic) {
1121       stochastic.setTemperature(targetTemperature);
1122     }
1123 
1124     // Set the step size.
1125     integrator.setTimeStep(dt);
1126 
1127     if (!initialized) {
1128       // Initialize from a restart file.
1129       if (loadRestart) {
1130         Crystal crystal = molecularAssembly[0].getCrystal();
1131         double[] x = state.x();
1132         double[] v = state.v();
1133         double[] a = state.a();
1134         double[] aPrevious = state.aPrevious();
1135         if (!dynFilter.readDYN(restartFile, crystal, x, v, a, aPrevious)) {
1136           String message = " Could not load the restart file - dynamics terminated.";
1137           logger.log(Level.WARNING, message);
1138           done = true;
1139           throw new IllegalStateException(message);
1140         } else {
1141           molecularAssembly[0].setCrystal(crystal);
1142         }
1143       } else {
1144         // Initialize using current atomic coordinates.
1145         potential.getCoordinates(state.x());
1146         // Initialize atomic velocities from a Maxwell-Boltzmann distribution or set to 0.
1147         if (initVelocities) {
1148           thermostat.maxwell(targetTemperature);
1149         } else {
1150           fill(state.v(), 0.0);
1151         }
1152       }
1153     } else {
1154       // If MD has already been run (i.e., Annealing or RepEx), then initialize velocities if
1155       // requested.
1156       if (initVelocities) {
1157         thermostat.maxwell(targetTemperature);
1158       }
1159     }
1160   }
1161 
1162   /** Initializes energy fields, esp. potential energy. */
1163   private void initializeEnergies() {
1164     // Compute the current potential energy.
1165     double[] x = state.x();
1166     double[] gradient = state.gradient();
1167 
1168     boolean propagateLambda = true;
1169     if (potential instanceof OrthogonalSpaceTempering orthogonalSpaceTempering) {
1170       propagateLambda = orthogonalSpaceTempering.getPropagateLambda();
1171       orthogonalSpaceTempering.setPropagateLambda(false);
1172     }
1173 
1174     if (esvSystem != null && potential instanceof OpenMMEnergy) {
1175       state.setPotentialEnergy(((OpenMMEnergy) potential).energyAndGradientFFX(x, gradient));
1176     } else {
1177       state.setPotentialEnergy(potential.energyAndGradient(x, gradient));
1178     }
1179 
1180     if (potential instanceof OrthogonalSpaceTempering orthogonalSpaceTempering) {
1181       orthogonalSpaceTempering.setPropagateLambda(propagateLambda);
1182     }
1183 
1184     // Initialize current and previous accelerations.
1185     if (!loadRestart || initialized || integrator instanceof Respa) {
1186       // For the Respa integrator, initial accelerations are from the slowly varying forces.
1187       if (integrator instanceof Respa) {
1188         potential.setEnergyTermState(Potential.STATE.SLOW);
1189         potential.energyAndGradient(x, gradient);
1190       }
1191 
1192       int numberOfVariables = state.getNumberOfVariables();
1193       double[] a = state.a();
1194       double[] mass = state.getMass();
1195       for (int i = 0; i < numberOfVariables; i++) {
1196         a[i] = -KCAL_TO_GRAM_ANG2_PER_PS2 * gradient[i] / mass[i];
1197       }
1198       state.copyAccelerationsToPrevious();
1199     }
1200 
1201     // Compute the current kinetic energy.
1202     thermostat.computeKineticEnergy();
1203 
1204     // Store the initial state.
1205     initialState = new UnmodifiableState(state);
1206 
1207     // Reset the statistics.
1208     temperatureStats.reset();
1209     potentialEnergyStats.reset();
1210     kineticEnergyStats.reset();
1211     totalEnergyStats.reset();
1212   }
1213 
1214   /** Pre-run operations (mostly logging) that require knowledge of system energy. */
1215   void postInitEnergies() {
1216     initialized = true;
1217 
1218     logger.log(basicLogging,
1219         format("\n  %8s %12s %12s %12s %8s %8s", "Time", "Kinetic", "Potential", "Total", "Temp", "CPU"));
1220     logger.log(basicLogging,
1221         format("  %8s %12s %12s %12s %8s %8s", "psec", "kcal/mol", "kcal/mol", "kcal/mol", "K", "sec"));
1222     logger.log(basicLogging, format("  %8s %12.4f %12.4f %12.4f %8.2f", "", state.getKineticEnergy(),
1223         state.getPotentialEnergy(), state.getTotalEnergy(), state.getTemperature()));
1224 
1225     // Store the initialized state.
1226     storeState();
1227   }
1228 
1229   /** Post-run cleanup operations. */
1230   void postRun() {
1231     // Add the potential energy of the slow degrees of freedom.
1232     if (integrator instanceof Respa) {
1233       potential.setEnergyTermState(Potential.STATE.BOTH);
1234     }
1235 
1236     // Log normal completion.
1237     if (!terminate) {
1238       logger.log(basicLogging, format(" Completed %8d time steps\n", nSteps));
1239     }
1240 
1241     // Reset the done and terminate flags.
1242     done = true;
1243     terminate = false;
1244   }
1245 
1246   /**
1247    * Append a snapshot to the trajectory file.
1248    *
1249    * @param extraLines Strings of meta-data to include.
1250    */
1251   protected void appendSnapshot(String[] extraLines) {
1252     // Loop over all molecular assemblies.
1253     for (MolecularAssembly assembly : molecularAssembly) {
1254       File archiveFile = assembly.getArchiveFile();
1255       ForceField forceField = assembly.getForceField();
1256       CompositeConfiguration properties = assembly.getProperties();
1257 
1258       // Save as an ARC file.
1259       if (archiveFile != null && !saveSnapshotAsPDB) {
1260         String aiName = FileUtils.relativePathTo(archiveFile).toString();
1261         if (esvSystem == null) {
1262           XYZFilter xyzFilter = new XYZFilter(archiveFile, assembly, forceField, properties);
1263           if (xyzFilter.writeFile(archiveFile, true, extraLines)) {
1264             logger.log(basicLogging, format(" Appended snapshot to:       %s", aiName));
1265           } else {
1266             logger.warning(          format(" Appending snapshot failed:  %s", aiName));
1267           }
1268         } else {
1269           XPHFilter xphFilter = new XPHFilter(archiveFile, assembly, forceField, properties,
1270               esvSystem);
1271           if (xphFilter.writeFile(archiveFile, true, extraLines)) {
1272             logger.log(basicLogging, format(" Appended to XPH archive %s", aiName));
1273           } else {
1274             logger.warning(format(" Appending to XPH archive %s failed.", aiName));
1275           }
1276         }
1277       } else if (saveSnapshotAsPDB) {
1278         File file = assembly.getFile();
1279         String extName = FilenameUtils.getExtension(file.getName());
1280         File pdbFile;
1281         if (extName.toLowerCase().startsWith("pdb")) {
1282           pdbFile = file;
1283         } else {
1284           String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
1285           pdbFile = new File(filename + ".pdb");
1286         }
1287         String aiName = FileUtils.relativePathTo(pdbFile).toString();
1288         PDBFilter pdbFilter = new PDBFilter(pdbFile, assembly, forceField, properties);
1289         if (pdbFilter.writeFile(pdbFile, true, extraLines)) {
1290           logger.log(basicLogging, format(" Appended to PDB file %s", aiName));
1291         } else {
1292           logger.warning(format(" Appending to PDB file to %s failed.", aiName));
1293         }
1294       }
1295     }
1296   }
1297 
1298   /**
1299    * Checks if thermodynamics must be logged. If logged, current time is returned, else the time
1300    * passed in is returned.
1301    *
1302    * @param step Time step to possibly log thermodynamics for.
1303    * @param time Clock time (in nsec) thermodynamics was last logged.
1304    * @return Either current time (if logged), else the time variable passed in.
1305    */
1306   protected long logThermoForTime(long step, long time) {
1307     if (step % logFrequency == 0) {
1308       return logThermodynamics(time);
1309     } else {
1310       return time;
1311     }
1312   }
1313 
1314   /**
1315    * Log thermodynamics to the screen.
1316    *
1317    * @param time Clock time (in nsec) since this was last called.
1318    * @return Clock time at the end of this call.
1319    */
1320   private long logThermodynamics(long time) {
1321     time = System.nanoTime() - time;
1322     logger.log(basicLogging,
1323         format(" %7.3e %12.4f %12.4f %12.4f %8.2f %8.3f", totalSimTime, state.getKineticEnergy(),
1324             state.getPotentialEnergy(), state.getTotalEnergy(), state.getTemperature(),
1325             time * NS2SEC));
1326     return System.nanoTime();
1327   }
1328 
1329   /** Main loop of the run method. */
1330   private void mainLoop() {
1331     // Integrate Newton's equations of motion for the requested number of steps,
1332     // unless early termination is requested.
1333     long time = System.nanoTime();
1334     for (long step = 1; step <= nSteps; step++) {
1335       if (step > 1) {
1336         List<Constraint> constraints = potential.getConstraints();
1337         // TODO: Replace magic numbers with named constants.
1338         long constraintFails = constraints.stream()
1339             .filter((Constraint c) -> !c.constraintSatisfied(state.x(), state.v(), 1E-7, 1E-7))
1340             .count();
1341         if (constraintFails > 0) {
1342           logger.info(format(" %d constraint failures in step %d", constraintFails, step));
1343         }
1344       }
1345 
1346       // Do the half-step thermostat operation.
1347       thermostat.halfStep(dt);
1348 
1349       // Do the half-step integration operation.
1350       integrator.preForce(potential);
1351       if (esvSystem != null) {
1352         esvIntegrator.preForce(potential);
1353         //preForce processes theta values after
1354         esvSystem.preForce();
1355       }
1356 
1357       // Compute the potential energy and gradients.
1358       if (esvSystem != null && potential instanceof OpenMMEnergy) {
1359         state.setPotentialEnergy(((OpenMMEnergy) potential).energyAndGradientFFX(state.x(), state.gradient()));
1360       } else {
1361         state.setPotentialEnergy(potential.energyAndGradient(state.x(), state.gradient()));
1362       }
1363 
1364       // Add the potential energy of the slow degrees of freedom.
1365       if (integrator instanceof Respa r) {
1366         double potentialEnergy = state.getPotentialEnergy();
1367         state.setPotentialEnergy(potentialEnergy + r.getHalfStepEnergy());
1368       }
1369 
1370       // Do the full-step integration operation.
1371       integrator.postForce(state.gradient());
1372       if (esvSystem != null) {
1373         double[] dEdL = esvSystem.postForce();
1374         esvIntegrator.postForce(dEdL);
1375       }
1376 
1377       // Compute the full-step kinetic energy.
1378       thermostat.computeKineticEnergy();
1379 
1380       // Do the full-step thermostat operation.
1381       thermostat.fullStep(dt);
1382 
1383       // Recompute the kinetic energy after the full-step thermostat operation.
1384       thermostat.computeKineticEnergy();
1385       if (esvSystem != null) {
1386         //Adiabatic thermostat does nothing at half step.
1387         esvThermostat.computeKineticEnergy();
1388       }
1389 
1390       // Remove center of mass motion every ~100 steps.
1391       int removeCOMMotionFrequency = 100;
1392       if (thermostat.getRemoveCenterOfMassMotion() && step % removeCOMMotionFrequency == 0) {
1393         thermostat.centerOfMassMotion(true, false);
1394       }
1395 
1396       // Collect current kinetic energy, temperature, and total energy.
1397       if (esvSystem != null) {
1398         double kineticEnergy = thermostat.getKineticEnergy();
1399         double esvKineticEnergy = esvThermostat.getKineticEnergy();
1400         state.setKineticEnergy(kineticEnergy + esvKineticEnergy);
1401       }
1402 
1403       // Collect running statistics.
1404       temperatureStats.addValue(state.getTemperature());
1405       potentialEnergyStats.addValue(state.getPotentialEnergy());
1406       kineticEnergyStats.addValue(state.getKineticEnergy());
1407       totalEnergyStats.addValue(state.getTotalEnergy());
1408 
1409       // Update atomic velocity, acceleration and previous acceleration.
1410       potential.setVelocity(state.v());
1411       potential.setAcceleration(state.a());
1412       potential.setPreviousAcceleration(state.aPrevious());
1413 
1414       // Log the current state every printFrequency steps.
1415       totalSimTime += dt;
1416       time = logThermoForTime(step, time);
1417       if (step % printEsvFrequency == 0 && esvSystem != null) {
1418         logger.log(basicLogging, format(" %s", esvSystem.getLambdaList()));
1419       }
1420 
1421       if (automaticWriteouts) {
1422         writeFilesForStep(step, true, true);
1423       }
1424 
1425       // Notify the algorithmListeners.
1426       if (algorithmListener != null && step % logFrequency == 0) {
1427         for (MolecularAssembly assembly : molecularAssembly) {
1428           algorithmListener.algorithmUpdate(assembly);
1429         }
1430       }
1431 
1432       // Check for a termination request.
1433       if (terminate) {
1434         logger.info(format("\n Terminating after %8d time steps\n", step));
1435         break;
1436       }
1437     }
1438   }
1439 
1440   /**
1441    * Perform a sanity check on a frequency to ensure it's not longer than total runtime. Currently,
1442    * the setter parameter is ignored due to issues with our test suite.
1443    *
1444    * @param describe Description of the frequency.
1445    * @param frequency Frequency in time steps.
1446    */
1447   private void checkFrequency(String describe, int frequency) {
1448     if (frequency > nSteps) {
1449       logger.fine(
1450           format(" Specified %s frequency of %d is greater than the number of steps %d", describe,
1451               frequency, nSteps));
1452     }
1453   }
1454 
1455   /**
1456    * Set the coordinates.
1457    *
1458    * @param coords The coordinates to copy into MD coordinates array.
1459    */
1460   public void setCoordinates(double[] coords) {
1461     state.setCoordinates(coords);
1462   }
1463 
1464   /**
1465    * Get the coordinates.
1466    *
1467    * @return A copy of the current coordinates are returned.
1468    */
1469   public double[] getCoordinates() {
1470     return state.getCoordinatesCopy();
1471   }
1472 
1473   /**
1474    * Store the current state of the molecular dynamics simulation in a MDState record.
1475    */
1476   public void storeState() {
1477     storedState = state.getUnmodifiableState();
1478   }
1479 
1480   /**
1481    * Revert the state of the MolecularDynamics instance to the stored MDState.
1482    *
1483    * @throws Exception is thrown if the stored state is null.
1484    */
1485   public void revertState() throws Exception {
1486     if (storedState == null) {
1487       throw new Exception();
1488     }
1489     revertState(storedState);
1490   }
1491 
1492   /**
1493    * Revert the state of the MolecularDynamics instance to the provided MDState.
1494    *
1495    * @param state The MDState to revert to.
1496    */
1497   private void revertState(UnmodifiableState state) {
1498     this.state.revertState(state);
1499     potential.setVelocity(state.v());
1500     potential.setAcceleration(state.a());
1501     potential.setPreviousAcceleration(state.aPrevious());
1502     // ToDo -- move these methods into the Potential interface.
1503     int numberOfVariables = this.state.getNumberOfVariables();
1504     Atom[] atoms = molecularAssembly[0].getActiveAtomArray();
1505     if (atoms.length * 3 == numberOfVariables) {
1506       int index = 0;
1507       double[] x = state.x();
1508       double[] gradient = state.gradient();
1509       for (Atom atom : atoms) {
1510         atom.moveTo(x[index], x[index + 1], x[index + 2]);
1511         atom.setXYZGradient(gradient[index], gradient[index + 1], gradient[index + 2]);
1512         index += 3;
1513       }
1514     }
1515   }
1516 
1517 }