AMS.cpp

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** SPDX-License-Identifier: BSD-3-Clause
**
** Copyright (c) 2010--present, eOn Development Team
** All rights reserved.
**
** Repo:
** https://github.com/TheochemUI/eOn
*/
 
#include "eon/potentials/AMS/AMS.h"
#include <cstddef>
#include <format>
#include <fstream>
#include <iostream>
#include <iterator>
#include <stdexcept>
#include <vector>
 
namespace bp = boost::process;
 
AMS::AMS(const Parameters &p)
    : Potential(PotType::AMS, p) {
  // Get the values from the configuration
  // All the parameter values convert to lowercase in generate_run
  this->engine = p.ams_options.engine;
  this->forcefield = p.ams_options.forcefield;
  this->model = p.ams_options.model;
  this->xc = p.ams_options.xc;
  this->resources = p.ams_options.resources;
  this->basis = p.ams_options.basis;
  this->engine_setup = generate_run(p);
  // Environment
  // TODO: Add more checks for how this can be set
  if (p.ams_options.env.amshome.empty() &&
      p.ams_options.env.scm_tmpdir.empty() &&
      p.ams_options.env.scmlicense.empty() &&
      p.ams_options.env.scm_pythondir.empty() &&
      p.ams_options.env.amsbin.empty() &&
      p.ams_options.env.amsresources.empty()) {
    nativenv = boost::this_process::environment();
  } else {
    nativenv = boost::this_process::environment();
    // Some of these can derive from the others
    nativenv["AMSHOME"] = p.ams_options.env.amshome;
    nativenv["SCM_TMPDIR"] = p.ams_options.env.scm_tmpdir;
    nativenv["SCMLICENSE"] = p.ams_options.env.scmlicense;
    nativenv["SCM_PYTHONDIR"] = p.ams_options.env.scm_pythondir;
    nativenv["AMSBIN"] = p.ams_options.env.amsbin;
    nativenv["AMSRESOURCES"] = p.ams_options.env.amsresources;
    nativenv["PATH"] += p.ams_options.env.amsbin;
  }
  // Do not pass "" in the config files
  // std::cout<<nativenv["PATH"].to_string()<<std::endl;
  amsevals = 0;
  // TODO: Optimize and reuse existing files Currently each Matter will
  // recreate the folders It should instead figure out if results exist and
  // use them
  this->first_run = true;
  // Determine if the engine supports restarts
  if (engine == "MOPAC") {
    this->can_restart = false;
    this->cjob = "amsResults";
    this->pjob = "amsResults";
  } else {
    this->can_restart = true;
    this->cjob = "firstRun";
    this->pjob = "secondRun";
  }
  return;
}
 
void AMS::cleanMemory(void) { return; }
 
AMS::~AMS() { cleanMemory(); }
 
namespace {
 
// The driver script, written afresh before every AMS invocation.
constexpr const char *kRunScript = "run_AMS.sh";
 
const char *elementArray[] = {
    "Unknown", "H",  "He", "Li", "Be", "B",  "C",  "N",  "O",  "F",  "Ne", "Na",
    "Mg",      "Al", "Si", "P",  "S",  "Cl", "Ar", "K",  "Ca", "Sc", "Ti", "V",
    "Cr",      "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br",
    "Kr",      "Rb", "Sr", "Y",  "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag",
    "Cd",      "In", "Sn", "Sb", "Te", "I",  "Xe", "Cs", "Ba", "La", "Ce", "Pr",
    "Nd",      "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu",
    "Hf",      "Ta", "W",  "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi",
    "Po",      "At", "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U",  NULL};
 
// guess the atom type from the atomic mass,
std::string mass2atom(double atomicmass) {
  return elementArray[int(atomicmass + .5)];
}
 
int symbol2atomicNumber(char const *symbol) {
  int i = 0;
 
  while (elementArray[i] != NULL) {
    if (strcmp(symbol, elementArray[i]) == 0) {
      return i;
    }
    i++;
  }
  // invalid symbol
  return -1;
}
 
char const *atomicNumber2symbol(int n) {
  // The trailing NULL terminates the table, so it bounds the valid range.
  if (n < 0 || static_cast<std::size_t>(n) + 1 >= std::size(elementArray)) {
    throw std::runtime_error(
        std::format("AMS knows no element symbol for atomic number {}", n));
  }
  return elementArray[n];
}
} // namespace
 
void AMS::runAMS() {
  boost::asio::io_context amsRun;
  std::future<std::string> run_out_future, run_err_future;
  std::string runout, runerr;
  if (chmod(kRunScript, S_IRWXU) != 0) {
    throw std::runtime_error(
        std::format("Could not make {} executable", kRunScript));
  }
  nativenv["AMS_JOBNAME"] = cjob;
  // assert(validate_order() == true);         // TODO: Debug only
  bp::child c(std::string(kRunScript),
              nativenv,                     // set the input
              bp::std_in.close(),           // no input
              bp::std_out > run_out_future, // STDOUT
              bp::std_err > run_err_future, // STDERR
              amsRun);
  amsRun.run();
  runout = run_out_future.get();
  runerr = run_err_future.get();
  if (runerr.find("ERROR") != std::string::npos) {
    bp::spawn("cat myrestart.in");
    bp::spawn("cat run_AMS.sh");
    throw std::runtime_error("\n AMS error while running");
  } else {
    this->amsevals = amsevals + 1;
    // std::cout << "Run completed normally" << std::endl;
    // std::cerr << "NORMAL TERMINATION\n";
  }
}
 
double AMS::extract_scalar_rkf(std::string key) {
  // The logic here handles the extraction asynchronously, as we lack
  // interest in processing this output line by line, and would prefer to
  // have it all in one place
  std::string execString;
  std::vector<std::string> execDat;
  boost::asio::io_context rkf;
  std::future<std::string> rkf_out_future, rkf_err_future;
  std::string rkfout, rkferr;
  double xval, x;
  std::vector<double> extracted;
  execString = std::format("dmpkf {}.results/{}.rkf AMSResults%{}", this->cjob,
                           this->engine_lower, key);
  // std::cout << execString << "\n";
  // Extract
  bp::child c(execString, nativenv,         // execute with the environment
              bp::std_in.close(),           // no input
              bp::std_out > rkf_out_future, // STDOUT
              bp::std_err > rkf_err_future, // STDERR
              rkf);
 
  rkf.run(); // this blocks until the end of the command
  // Populate the strings
  rkfout = rkf_out_future.get();
  rkferr = rkf_err_future.get();
  if (rkferr.find("ERROR") != std::string::npos) {
    throw std::runtime_error(std::format(
        "\n AMS error while extracting {}, got:\n {}", key, rkferr));
  } else {
    // std::cout << "Extracting " << key << std::endl;
  }
 
  execDat = absl::StrSplit(rkfout, '\n');
  // Is energy or some other scalar property
  // First 3 lines hold the headers
  // [0] = "AMSResults                       "
  // [1] = "Energy                           "
  // [2] = "         1         1         2"
  // [3] = "   0.135012547958282714E+000"
  // [4] = ""
  if (execDat.size() < 4) {
    throw std::runtime_error(
        std::format("\n AMS returned {} lines for {}, too few to hold a value",
                    execDat.size(), key));
  }
  if (absl::SimpleAtod(execDat[3], &x)) {
    xval = x * this->energyConversion;
    // std::cout << std::format(
    //     "\n Got {:.4f} Hartree from AMS and converted to {:.4f} eV\n", x,
    //     xval);
    return xval;
  } else {
    throw std::runtime_error(
        std::format("\n Expected {}, got {} instead", key, execDat[3]));
  }
  // Never reach here
  throw std::runtime_error("Generic AMS dmpkf scalar error \n");
}
 
std::vector<double> AMS::extract_cartesian_rkf(std::string key) {
  std::string execString;
  std::vector<std::string> execDat, innerDat;
  boost::asio::io_context rkf;
  std::future<std::string> rkf_out_future, rkf_err_future;
  std::string rkfout, rkferr;
  double felem, x;
  std::vector<double> extracted;
  execString = std::format("dmpkf {}.results/{}.rkf AMSResults%{}", this->cjob,
                           this->engine_lower, key);
  // std::cout << execString << "\n";
 
  // Extract
  bp::child c(execString, nativenv,         // execute with the environment
              bp::std_in.close(),           // no input
              bp::std_out > rkf_out_future, // STDOUT
              bp::std_err > rkf_err_future, // STDERR
              rkf);
 
  rkf.run(); // this blocks until the end of the command
  // Populate the strings
  rkfout = rkf_out_future.get();
  rkferr = rkf_err_future.get();
 
  if (rkferr.find("ERROR") != std::string::npos) {
    throw std::runtime_error(std::format(
        "\n AMS error while extracting {}, got:\n {}", key, rkferr));
  } else {
    // std::cout << "Extracting " << key << std::endl;
  }
 
  execDat = absl::StrSplit(rkfout, '\n');
  // Assume gradients or some other x y z property
  for (int i = 3; i < execDat.size(); i++) {
    // There exist one per atom, the number of which equals N
    // The first three lines hold the header files as before
    std::vector<std::string> strrow = absl::StrSplit(execDat[i], ' ');
    for (auto elem : strrow) {
      // This loop uses an auto variable since some of the elements of strrow
      // often evaluate to ' '
      // TODO: Optimize, perhaps with a regex
      // [0] = ""
      // [1] = ""
      // [2] = ""
      // [3] = "0.798885933949825835E-004"
      // [4] = ""
      // [5] = "-0.755400038198822616E-004"
      // [6] = ""
      // [7] = ""
      // [8] = "0.457005740252359742E-004"
      if (!elem.empty() && absl::SimpleAtod(elem, &x)) {
        felem = x * this->forceConversion;
        // std::cout << std::format(
        //     "\n Gradient element={:.4f} Hartree/Bohr from AMS\n Force "
        //     "element={:.4f} eV/Angstrom\n", x, felem);
        extracted.emplace_back(felem);
      }
    }
  }
  // // Debug
  // int counter = 0;
  // for (int a = 0; a < N * 3; a++) {
  //   std::cout << std::format("{:.25e} ", forces[a]);
  //   counter++;
  //   if (counter % 3 == 0) {
  //     std::cout << std::endl;
  //   }
  // }
  return extracted;
}
 
void AMS::updateCoord(long N, const double *R) {
  // The logic used here requires us to update the PREVIOUS job, since that
  // functions as the one from which the calculation restarts
  // Only the third line presents vague complexity
  // https://www.scm.com/doc/Scripting/Commandline_Tools/KF_command_line_utilities.html
  std::ofstream updCoord;
  std::string execString, coordDump, newCoord;
  std::vector<std::string> execDat;
  boost::asio::io_context coordio;
  std::future<std::string> err, rdump;
  std::vector<double> gradients;
  // Prep new run
  // Get the previous run's coordinates
  execString = std::format("dmpkf {}.results/ams.rkf Molecule%Coords", pjob);
  // std::cout << execString << "\n";
  // Store Coordinates
  // TODO: Simplify this, we only need the first few lines
  bp::child cprog(execString, nativenv, bp::std_in.close(), bp::std_out > rdump,
                  bp::std_err > err, coordio);
  coordio.run();
  execDat = absl::StrSplit(rdump.get(), '\n');
  // Get the first three lines
  int counter = 0;
  for (auto j : execDat) {
    if (counter >= 3) {
      break;
    } else {
      absl::StrAppend(&newCoord, j, "\n");
    }
    counter++;
  }
  // Put the rest of the coordinates
  counter = 0;
  for (int a = 0; a < N * 3; a++) {
    absl::StrAppend(&newCoord, (R[a] * lengthConversion), "  ");
    counter++;
    if (counter % 3 == 0) {
      absl::StrAppend(&newCoord, "\n");
    }
  }
  coordDump = "#!/bin/sh\n udmpkf ";
  absl::StrAppend(&coordDump, pjob, ".results/ams.rkf <<EOF\n", newCoord,
                  "EOF");
  // std::cout << coordDump;
  updCoord.open("updCoord.sh", std::ios::trunc);
  if (!updCoord) {
    throw std::runtime_error("Could not open updCoord.sh for writing");
  }
  updCoord << coordDump;
  updCoord.close();
  if (!updCoord) {
    throw std::runtime_error("Could not write the coordinates to updCoord.sh");
  }
  // bp::spawn("chmod +x updCoord.sh");
  if (chmod("updCoord.sh", S_IRWXU) != 0) {
    throw std::runtime_error("Could not make updCoord.sh executable");
  }
  bp::child cuprog("updCoord.sh", nativenv, bp::std_err > bp::null);
  cuprog.wait();
  return;
}
 
void AMS::switchjob() {
  std::string tmp;
  // std::cout << std::format("\nEntered Switch:\nCurrent:{}, Previous:{}\n",
  // cjob, pjob);
  tmp = this->cjob;
  this->cjob = this->pjob;
  this->pjob = tmp;
  // std::cout << std::format("\nSwitched\n Current:{}, Previous:{}\n", cjob,
  // pjob);
}
 
void AMS::write_restart() {
  std::string restart_formatter;
  if (can_restart) {
    restart_formatter = R"(
   EngineRestart {0}.results/{1}.rkf
   LoadSystem
    File {0}.results/ams.rkf
    Section Molecule
   End
  )";
  } else {
    restart_formatter = R"(
   LoadSystem
    File {0}.results/ams.rkf
    Section Molecule
   End
  )";
  }
  std::string restart_data = std::format(restart_formatter, pjob, engine_lower);
  restartFrom.open("myrestart.in", std::ios::trunc);
  if (!restartFrom) {
    throw std::runtime_error("Could not open myrestart.in for writing");
  }
  restartFrom << restart_data;
  restartFrom.close();
  if (!restartFrom) {
    throw std::runtime_error("Could not write the restart block to "
                             "myrestart.in");
  }
  return;
}
 
void AMS::copyForces(long N, const std::vector<double> &frc, double *F) {
  // A short AMS run yields fewer than 3N values; copying blindly would hand
  // the optimizer whatever lies past the end of the vector.
  if (frc.size() < static_cast<std::size_t>(3 * N)) {
    throw std::runtime_error(
        std::format("\n AMS returned {} gradient components, expected {}",
                    frc.size(), 3 * N));
  }
  for (long i = 0; i < 3 * N; i++) {
    F[i] = frc[i];
  }
}
 
void AMS::force(long N, const double *R, const int *atomicNrs, double *F,
                double *U, double *variance, const double *box) {
  variance = nullptr;
  if (not can_restart or first_run) {
    // std::cout << std::format("\nCAN_RESTART:{}  FIRST_RUN:{}\n", can_restart,
    // first_run); This holds true for all engines with no restart Also if an
    // engine supports being restarted, the first run needs this
    passToSystem(N, R, atomicNrs, box);
    runAMS();
    copyForces(N, extract_cartesian_rkf("Gradients"), F);
    *U = extract_scalar_rkf("Energy");
    // Update, will still not matter for those without restarts
    if (can_restart) {
      first_run = false;
      // We need the "previous job" to pre-populate
      // clang-format off
      // std::cout<<std::format("\nMoving {} to {} before switching during the first job\n", cjob, pjob);
      const auto copyOptions = std::filesystem::copy_options::overwrite_existing
                             | std::filesystem::copy_options::recursive
                             ;
      // clang-format on
      std::filesystem::copy(std::format("./{}.results/", cjob),
                            std::format("./{}.results/", pjob), copyOptions);
    }
    return;
  } else {
    // std::cout << std::format("\nCAN_RESTART:{}  FIRST_RUN:{}\n", can_restart,
    // first_run);
    smallSys(N, R, atomicNrs, box); // writes run_AMS.sh
    updateCoord(N, R);              // updates coordinates in previous job
    write_restart();                // writes restart file using previous job
    runAMS();
    copyForces(N, extract_cartesian_rkf("Gradients"), F);
    *U = extract_scalar_rkf("Energy");
    switchjob(); // toggles the jobs
    return;
  }
  // Never reach here
  throw std::runtime_error("Generic AMS force error \n");
}
 
void AMS::passToSystem(long N, const double *R, const int *atomicNrs,
                       const double *box)
// Creating the standard input file that the AMS driver reads
{
  std::ofstream out(kRunScript, std::ios::trunc);
  if (!out) {
    throw std::runtime_error(
        std::format("Could not open {} for writing", kRunScript));
  }
 
  out << "#!/bin/sh\n";
  out << std::format("export AMS_JOBNAME={}\n", cjob);
  out << "$AMSBIN/ams --delete-old-results <<eor\n";
  out << "Task SinglePoint\n";
  out << "System\n";
  out << " Atoms\n";
  for (long i = 0; i < N; i++) {
    out << std::format("  {}\t{:.19f}\t{:.19f}\t{:.19f}\n",
                       atomicNumber2symbol(atomicNrs[i]), R[i * 3 + 0],
                       R[i * 3 + 1], R[i * 3 + 2]);
  }
  out << " End\n";
  if (not model.empty() || not forcefield.empty()) {
    out << " Lattice\n";
    for (int i = 0; i < 3; i++) {
      out << std::format("  {:.19f}\t{:.19f}\t{:.19f}\n", box[i * 3 + 0],
                         box[i * 3 + 1], box[i * 3 + 2]);
    }
    out << " End\n";
  }
  out << "End\n";
  out << engine_setup;
  out << "Properties\n";
  out << " Gradients\n";
  out << "End\n";
  if (can_restart and not first_run) {
    out << "@include myrestart.in\n";
  }
  out << "eor";
 
  finishRunScript(out);
  return;
}
 
void AMS::smallSys(long N, const double *R, const int *atomicNrs,
                   const double *box)
// Creating the truncated input file that the AMS driver reads
{
  std::ofstream out(kRunScript, std::ios::trunc);
  if (!out) {
    throw std::runtime_error(
        std::format("Could not open {} for writing", kRunScript));
  }
 
  out << "#!/bin/sh\n";
  out << std::format("export AMS_JOBNAME={}\n", cjob);
  out << "$AMSBIN/ams --delete-old-results <<eor\n";
  out << "Task SinglePoint\n";
  out << engine_setup;
  out << "Properties\n";
  out << " Gradients\n";
  out << "End\n";
  out << "@include myrestart.in\n";
  out << "eor";
 
  finishRunScript(out);
  return;
}
 
void AMS::finishRunScript(std::ofstream &out) {
  out.close();
  if (!out) {
    throw std::runtime_error(
        std::format("Could not write the AMS input to {}", kRunScript));
  }
  if (chmod(kRunScript, S_IRWXU) != 0) {
    throw std::runtime_error(
        std::format("Could not make {} executable", kRunScript));
  }
}
 
std::string AMS::generate_run(const Parameters &p) {
  std::string engine_block; // Shadows the class variable
  // TODO: Use args everywhere, cleaner logic
  // Ensure capitals and existence
  engine.empty()
      ? throw std::runtime_error("AMS Engine is required \n")
      : std::transform(engine.begin(), engine.end(), engine.begin(), ::toupper);
  // engine functions uniquely, it serves as a filename, so we store
  // engine_lower as a lowercase version too
  engine_lower = engine;
  std::transform(engine.begin(), engine.end(), engine_lower.begin(), ::tolower);
  // Prepare the block
  if (engine == "MOPAC") {
    model.empty()
        ? throw std::runtime_error("MOPAC needs a model\n")
        : std::transform(model.begin(), model.end(), model.begin(), ::toupper);
    std::string engine_formatter = R"(
 Engine {}
  Model {}
 EndEngine
)";
    engine_block = std::format(engine_formatter, engine, model);
    return engine_block;
  } else if (engine == "ADF" || engine == "BAND") {
    basis.empty()
        ? throw std::runtime_error("ADF/BAND need a basis\n")
        : std::transform(basis.begin(), basis.end(), basis.begin(), ::toupper);
    xc.empty() ? throw std::runtime_error("ADF/BAND need a functional\n")
               : std::transform(xc.begin(), xc.end(), xc.begin(), ::toupper);
    std::string engine_formatter = R"(
   Engine {}
     Basis
       Type {}
     End
 
     XC
       {}
     End
   EndEngine
  )";
    engine_block = std::format(engine_formatter, engine, basis, xc);
    return engine_block;
  } else if (engine == "DFTB") {
    resources.empty() ? throw std::runtime_error("DFTB need resources\n")
                      : std::transform(resources.begin(), resources.end(),
                                       resources.begin(), ::toupper);
    std::string engine_formatter = R"(
   Engine {}
     ResourcesDir {}
   EndEngine
  )";
    engine_block = std::format(engine_formatter, engine, resources);
    return engine_block;
  } else if (engine == "reaxff") {
    forcefield.empty() ? throw std::runtime_error("REAXFF needs a forcefield\n")
                       : std::transform(forcefield.begin(), forcefield.end(),
                                        forcefield.begin(), ::toupper);
 
    std::string engine_formatter = R"(
   Engine {}
     ForceField {}
   EndEngine
  )";
    engine_block = std::format(engine_formatter, engine, forcefield);
    return engine_block;
  } else if (engine == "FORCEFIELD") {
    std::string engine_formatter = R"(
   Engine {}
   EndEngine
  )";
    engine_block = std::format(engine_formatter, engine);
    return engine_block;
  }
 
  // Never reach here
  throw std::runtime_error("Generic AMS engine error \n");
}
 
/*
**
** DEBUGGER
**
** The set of functions below, when activated, include asserts to ensure
*equality with AMS_IO
* These work.
*/
 
// clang-format off
 
// void AMS::recieveFromSystem(long N, double *F, double *U) {
 
//   FILE *in;
//   double junkF;
//   char junkChar[256];
//   double forceX;
//   double forceY;
//   double forceZ;
//   double index;
//   char line[256];
 
//   in = fopen("ams_output", "r");
 
//   while (fgets(line, sizeof(line), in)) {
 
//     if (strcmp(line, "      CALCULATION RESULTS\n") ==
//         0) { // Finding the Energy in the output file
 
//       fscanf(in, "%s %s %s %lf", junkChar, junkChar, junkChar, U);
//       *U = *U * 27.2114; // Energy in hartree to eV
//     }
//     if (strcmp(line, "  Index   Atom            d/dx            d/dy           "
//                      " d/dz\n") == 0) { // Finding the forces
//       for (int i = 0; i < N; i++) {
//         fscanf(in, "%lf %s %lf %lf %lf", &index, &junkChar, &forceX, &forceY,
//                &forceZ);
//         F[int(i) * 3 + 0] = -forceX;
//         F[int(i) * 3 + 1] = -forceY;
//         F[int(i) * 3 + 2] = -forceZ;
//         // AMS gives gradients, not forces, hence the change.
//       }
//     }
//   }
//   for (int i = 0; i < 3 * N; i++) {
//     F[i] = F[i] * 51.4220862; // Forces from hartree/bohr to eV/Angstrom
//   }
 
//   fclose(in);
//   return;
// }
 
// void AMS::force(long N, const double *R, const int *atomicNrs, double *F,
//                 double *U, const double *box, int nImages = 1) {
//   passToSystem(N, R, atomicNrs, box);
//   system("chmod +x run_AMS.sh");
//   system("./run_AMS.sh >> ams_output"); // Run a single point AMS calculation
//                                         // and write the results into ams_output
//   recieveFromSystem(N, F, U);
//   runAMS();
//   double energ = extract_scalar_rkf("Energy"); // Sets energy
//   // assert(fabs(*U-this->energy)<DBL_EPSILON); // Equality == doesn't work well
//   // for floats
//   assert(fabs(*U - energ) < 1e-5); // Equality == doesn't work well for floats
//   // *U = this->energy;
//   // this->forces.clear();        // TODO: Slow!
//   std::vector<double> frc = extract_cartesian_rkf("Gradients"); // Sets forces
//   double *ftest = frc.data();
//   for (int i = 0; i < forces.size(); i++) {
//     assert(fabs(F[i] - ftest[i]) < 1e-5);
//   }
//   *U = energ;
//   F = ftest;
//   return;
// }
 
// Ordering which works, scopes are correct
// void AMS::force(long N, const double *R, const int *atomicNrs, double *F,
//                 double *U, const double *box, int nImages = 1) {
//   // Somehow broken, even though this matches the prior configuration
//   passToSystem(N, R, atomicNrs, box);
//   runAMS();
//   std::vector<double> frc = extract_cartesian_rkf("Gradients");
//   double *ftest = frc.data();
//   for (int i = 0; i < N; i++) {
//     F[3 * i] = ftest[3 * i];
//     F[3 * i + 1] = ftest[3 * i + 1];
//     F[3 * i + 2] = ftest[3 * i + 2];
//   }
//   *U = extract_scalar_rkf("Energy");
//   return;
// }
 
// clang-format on
 
/*
** Debugging Toggles
** These functions validate the job ordering when added to
*runAMS()
*/
 
// TODO: Only in debug
// std::string AMS::readFile(std::filesystem::path path) {
//   // Kanged: https://stackoverflow.com/a/40903508/1895378
//   std::ifstream f(path, std::ios::in | std::ios::binary);
//   const auto sz = std::filesystem::file_size(path);
//   std::string stres(sz, '\0');
//   f.read(stres.data(), sz);
//   return stres;
// }
 
// bool AMS::validate_order() {
//   // Validate all inputs
//   if (not first_run) {
//     std::string rfile = readFile("run_AMS.sh");
//     std::string resfile = readFile("myrestart.in");
//     std::string updcoord = readFile("updCoord.sh");
//     // ifstream runfile("run_AMS.sh"), restartfile("myrestart.in"),
//     // updcoord("updCoord.sh"); istream_iterator<string> rfiter(runfile),
//     // refiter(restartfile), uciter(updcoord), eof; vector<string>
//     // rfstore(rfiter, eof), refstore(refiter, eof), ucstore(uciter, eof);
//     // (?<=AMS_JOBNAME=).*$
//     // (?<=udmpkf ).*(?=\.results)
//     // (?<=File ).*(?=\.results)
//     // The logic here enforces that the current job restarts from the
//     previous one, so
//     // the coordinates of the previous job update, the restart
//     references
//     // the previous job, and then finally the current job executes with cjob
//     assert(absl::StrContains(rfile, cjob));
//     assert(absl::StrContains(resfile, pjob));
//     assert(absl::StrContains(updcoord, pjob));
//   }
//   return true;
// }