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 "eon/Parameters.h"
#include <algorithm>
#include <cctype>
#include <cstddef>
#include <cstdint>
#include <format>
#include <fstream>
#include <iterator>
#include <ranges>
#include <readcon-core.hpp>
#include <stdexcept>
#include <vector>
 
namespace bp = boost::process;
 
AMS::AMS(const eonc::Parameters &p)
    : eonc::Potential(eonc::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
  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";
 
std::string symbol_for_z(int n) {
  if (n <= 0) {
    throw std::runtime_error(
        std::format("AMS knows no element symbol for atomic number {}", n));
  }
  return readcon::z_to_symbol(static_cast<uint64_t>(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;
  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;
  }
}
 
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);
  // 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));
  }
 
  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;
    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);
 
  // 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));
  }
 
  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;
        extracted.emplace_back(felem);
      }
    }
  }
  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);
  // 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");
  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;
  tmp = this->cjob;
  this->cjob = this->pjob;
  this->pjob = tmp;
}
 
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) {
    // 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
      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 {
    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",
                       symbol_for_z(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 eonc::Parameters &p) {
  std::string engine_block; // Shadows the class variable
  // TODO: Use args everywhere, cleaner logic
  // Ensure capitals and existence
  if (engine.empty()) {
    throw std::runtime_error("AMS Engine is required \n");
  }
  std::ranges::transform(engine, engine.begin(), [](unsigned char c) {
    return static_cast<char>(std::toupper(c));
  });
  // engine functions uniquely, it serves as a filename, so we store
  // engine_lower as a lowercase version too
  engine_lower = engine;
  std::ranges::transform(engine, engine_lower.begin(), [](unsigned char c) {
    return static_cast<char>(std::tolower(c));
  });
  // Prepare the block
  if (engine == "MOPAC") {
    if (model.empty()) {
      throw std::runtime_error("MOPAC needs a model\n");
    }
    std::ranges::transform(model, model.begin(), [](unsigned char c) {
      return static_cast<char>(std::toupper(c));
    });
    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") {
    if (basis.empty()) {
      throw std::runtime_error("ADF/BAND need a basis\n");
    }
    std::ranges::transform(basis, basis.begin(), [](unsigned char c) {
      return static_cast<char>(std::toupper(c));
    });
    if (xc.empty()) {
      throw std::runtime_error("ADF/BAND need a functional\n");
    }
    std::ranges::transform(xc, xc.begin(), [](unsigned char c) {
      return static_cast<char>(std::toupper(c));
    });
    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") {
    if (resources.empty()) {
      throw std::runtime_error("DFTB need resources\n");
    }
    std::ranges::transform(resources, resources.begin(), [](unsigned char c) {
      return static_cast<char>(std::toupper(c));
    });
    std::string engine_formatter = R"(
   Engine {}
     ResourcesDir {}
   EndEngine
  )";
    engine_block = std::format(engine_formatter, engine, resources);
    return engine_block;
  } else if (engine == "reaxff") {
    if (forcefield.empty()) {
      throw std::runtime_error("REAXFF needs a forcefield\n");
    }
    std::ranges::transform(forcefield, forcefield.begin(), [](unsigned char c) {
      return static_cast<char>(std::toupper(c));
    });
 
    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");
}