NudgedElasticBandJob.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/NudgedElasticBandJob.h"
#include "eon/ConjugateGradients.h"
#include "eon/EonLogger.h"
#include "eon/NEBInitialPaths.hpp"
#include "eon/NEBSplineExtrema.h"
#include "eon/Potential.h"
 
#include <filesystem>
#include <format>
#include <fstream>
#include <stdexcept>
 
std::vector<std::string> NudgedElasticBandJob::run() {
  NudgedElasticBand::NEBStatus status;
  size_t f1;
 
  std::string reactantFilename = eonc::helpers::getRelevantFile("reactant.con");
  std::string productFilename = eonc::helpers::getRelevantFile("product.con");
 
  // TS interpolation: only when initializer is FILE (explicit path loading).
  // Previously, any ts.con in CWD would silently override SIDPP/IDPP paths.
  bool tsInterpolate = false;
  std::shared_ptr<Matter> transitionState = nullptr;
  if (params.neb_options.initialization.method == NEBInit::FILE &&
      std::filesystem::exists("ts.con")) {
    std::string transitionStateFilename =
        eonc::helpers::getRelevantFile("ts.con");
    transitionState = std::make_shared<Matter>(pot, params);
    tsInterpolate = true;
    if (!eonc::io::io_ok(
            transitionState->con2matter(transitionStateFilename))) {
      EONC_LOG_CRITICAL("Failed to load {}", transitionStateFilename);
      throw std::runtime_error("failed to load " + transitionStateFilename);
    }
  }
 
  auto initial = std::make_shared<Matter>(pot, params);
  auto final_state = std::make_shared<Matter>(pot, params);
 
  // When initializer=file and initial_path_in is set, endpoints are the first
  // and last frames of that list — do not require reactant.con / product.con
  // (issue #278; buggy endpoint files were still loaded and could crash).
  const bool usePathEndpoints =
      params.neb_options.initialization.method == NEBInit::FILE &&
      !params.neb_options.initialization.input_path.empty();
 
  if (usePathEndpoints) {
    const auto file_paths = eonc::helpers::neb_paths::readFilePaths(
        params.neb_options.initialization.input_path);
    if (file_paths.size() < 2) {
      throw std::runtime_error(
          "NEB initial_path_in must list at least two frames "
          "(reactant and product endpoints)");
    }
    if (!eonc::io::io_ok(initial->con2matter(file_paths.front().string()))) {
      EONC_LOG_CRITICAL("Failed to load NEB path reactant frame {}",
                        file_paths.front().string());
      throw std::runtime_error("failed to load NEB path reactant frame");
    }
    if (!eonc::io::io_ok(final_state->con2matter(file_paths.back().string()))) {
      EONC_LOG_CRITICAL("Failed to load NEB path product frame {}",
                        file_paths.back().string());
      throw std::runtime_error("failed to load NEB path product frame");
    }
    QUILL_LOG_INFO(m_log,
                   "NEB endpoints from initial path list (skipped {} / {})",
                   reactantFilename, productFilename);
  } else {
    if (!eonc::io::io_ok(initial->con2matter(reactantFilename))) {
      EONC_LOG_CRITICAL("Failed to load {}", reactantFilename);
      throw std::runtime_error("failed to load reactant.con");
    }
    if (!eonc::io::io_ok(final_state->con2matter(productFilename))) {
      EONC_LOG_CRITICAL("Failed to load {}", productFilename);
      throw std::runtime_error("failed to load product.con");
    }
  }
 
  // Endpoint minimization logic:
  // - If params.neb_options.endpoints.minimize is false: never minimize
  // endpoints.
  // - If params.neb_options.endpoints.minimize is true and
  // params.neb_options.initialization.input_path is empty: minimize endpoints.
  // - If params.nebMinimEP is true and params.nebIpath is NOT empty:
  //     -> minimize endpoints only if params.nebMinimEPIpath is true.
  // Log what decision was made so users can see behavior.
  bool shouldMinimizeEndpoints = false;
  if (!params.neb_options.endpoints.minimize) {
    QUILL_LOG_DEBUG(m_log,
                    "minimize_endpoints == false: not minimizing endpoints.");
    shouldMinimizeEndpoints = false;
  } else {
    if (params.neb_options.initialization.input_path.empty()) {
      QUILL_LOG_DEBUG(m_log, "minimize_endpoints == true and nebIpath "
                             "empty: minimizing endpoints.");
      shouldMinimizeEndpoints = true;
    } else {
      // nebIpath provided: only minimize if neb_options.endpoints.use_path_file
      // explicitly allowed.
      if (params.neb_options.endpoints.use_path_file) {
        QUILL_LOG_DEBUG(
            m_log,
            "minimize_endpoints == true and nebIpath provided, but "
            "minimize_endpoints_for_ipath == true: minimizing endpoints.");
        shouldMinimizeEndpoints = true;
      } else {
        QUILL_LOG_DEBUG(
            m_log,
            "minimize_endpoints == true but nebIpath provided and "
            "minimize_endpoints_for_ipath == false: not minimizing endpoints.");
        shouldMinimizeEndpoints = false;
      }
    }
  }
 
  if (shouldMinimizeEndpoints) {
    QUILL_LOG_DEBUG(m_log, "Minimizing reactant");
    // TODO(rg): Maybe when we have even more parameters, false can be set by
    // the user too..
    initial->relax(false, params.debug_options.write_movies,
                   params.main_options.checkpoint, "react_neb", "react_neb");
    // TODO(rg): How do we report the total E/F now? Currently this is just the
    // total total, people might want "per-stage" totals (but they can also get
    // them from the log.)
    QUILL_LOG_DEBUG(m_log, "Minimized reactant in ");
    QUILL_LOG_DEBUG(m_log, "Minimizing product");
    final_state->relax(false, params.debug_options.write_movies,
                       params.main_options.checkpoint, "prod_neb", "prod_neb");
  }
 
  auto neb =
      std::make_unique<NudgedElasticBand>(initial, final_state, params, pot);
 
  if (tsInterpolate) {
    AtomMatrix reactantToTS = transitionState->pbc(
        transitionState->getPositions() - initial->getPositions());
    AtomMatrix TSToProduct = transitionState->pbc(
        final_state->getPositions() - transitionState->getPositions());
    for (int image = 1; image <= neb->numImages; image++) {
      int mid = neb->numImages / 2 + 1;
      if (image < mid) {
        double frac = static_cast<double>(image) / static_cast<double>(mid);
        neb->path[image]->setPositions(initial->getPositions() +
                                       frac * reactantToTS);
      } else if (image > mid) {
        double frac = static_cast<double>(image - mid) /
                      static_cast<double>(neb->numImages - mid + 1);
        neb->path[image]->setPositions(transitionState->getPositions() +
                                       frac * TSToProduct);
      } else if (image == mid) {
        neb->path[image]->setPositions(transitionState->getPositions());
      }
    }
  }
 
  f1 = PotRegistry::get().total_force_calls();
  status = neb->compute();
  fCallsNEB += PotRegistry::get().total_force_calls() - f1;
 
  if (status == NudgedElasticBand::NEBStatus::GOOD) {
    neb->printImageData();
    neb->findExtrema();
  }
 
  printEndState(status);
  saveData(status, neb.get());
 
  return returnFiles;
}
 
void NudgedElasticBandJob::saveData(NudgedElasticBand::NEBStatus status,
                                    NudgedElasticBand *neb) {
  std::string resultsFilename("results.dat");
  returnFiles.push_back(resultsFilename);
 
  {
    std::ofstream out(resultsFilename, std::ios::binary);
    if (!out) {
      QUILL_LOG_ERROR(m_log, "Failed to open {} for writing", resultsFilename);
      return;
    }
 
    out << std::format("{} termination_reason\n", static_cast<int>(status));
    out << std::format("{} termination_reason_text\n",
                       magic_enum::enum_name(status));
    out << std::format(
        "{} potential_type\n",
        magic_enum::enum_name<PotType>(params.potential_options.potential));
    out << std::format("{} total_force_calls\n",
                       PotRegistry::get().total_force_calls());
    out << std::format("{} force_calls_neb\n", fCallsNEB);
    out << std::format("{:f} energy_reference\n",
                       neb->path[0]->getPotentialEnergy());
    out << std::format("{} number_of_images\n", neb->numImages);
 
    for (long i = 0; i <= neb->numImages + 1; i++) {
      out << std::format("{:f} image{}_energy\n",
                         neb->path[i]->getPotentialEnergy() -
                             neb->path[0]->getPotentialEnergy(),
                         i);
      out << std::format("{:f} image{}_force\n",
                         neb->path[i]->getForces().norm(), i);
      double proj_norm = (i >= 1 && i <= neb->numImages)
                             ? neb->projectedForce[i]->norm()
                             : 0.0;
      out << std::format("{:f} image{}_projected_force\n", proj_norm, i);
    }
 
    long safeNumExtrema = std::min(
        neb->numExtrema, static_cast<long>(neb->extremumPosition.size()));
    out << std::format("{} number_of_extrema\n", safeNumExtrema);
    for (long i = 0; i < safeNumExtrema; i++) {
      out << std::format("{:f} extremum{}_position\n", neb->extremumPosition[i],
                         i);
      out << std::format("{:f} extremum{}_energy\n", neb->extremumEnergy[i], i);
    }
  }
 
  // Save the Full NEB Path
  std::string nebFilename("neb.con");
  returnFiles.push_back(nebFilename);
  if (!eonc::io::io_ok(eonc::neb::writePathCon(
          neb->path, neb->tangent, neb->eigenmode_solvers, neb->numImages,
          params.debug_options.estimate_neb_eigenvalues, nebFilename))) {
    QUILL_LOG_ERROR(m_log, "Failed to write {}", nebFilename);
  }
 
  // Save Discrete Saddle Point (Highest Energy Image)
  std::string spFilename("sp.con");
  if (!eonc::io::io_ok(
          neb->path[neb->maxEnergyImage]->matter2con(spFilename))) {
    QUILL_LOG_ERROR(m_log, "Failed to write {}", spFilename);
  }
  returnFiles.push_back(spFilename);
 
  // Setup Dimer Configurations for Each Spline Peak
  if (params.neb_options.mmf_peaks.enabled && neb->numExtrema > 0) {
    int peakCount = 0;
    for (long i = 0; i < neb->numExtrema; i++) {
      // Filter 1: Only look at maxima (negative curvature)
      // Filter 2: Energy threshold (e.g., peak must be > 0.05 eV above
      // reactant)
      double relativeEnergy =
          neb->extremumEnergy[i] - neb->path[0]->getPotentialEnergy();
 
      if (neb->extremumCurvature[i] < 0 &&
          relativeEnergy > params.neb_options.mmf_peaks.tolerance) {
        double posFraction = neb->extremumPosition[i];
        int leftIdx = static_cast<int>(std::floor(posFraction));
        double f = posFraction - leftIdx;
 
        if (leftIdx < 0 || leftIdx >= neb->numImages + 1)
          continue;
 
        // 1. Write Interpolated Position (.con)
        Matter peakPos = eonc::helpers::neb_paths::interpolateImage(
            *neb->path[leftIdx], *neb->path[leftIdx + 1], f);
        std::string peakPosFile = std::format("peak{:02d}_pos.con", peakCount);
        if (!eonc::io::io_ok(peakPos.matter2con(peakPosFile))) {
          QUILL_LOG_ERROR(m_log, "Failed to write {}", peakPosFile);
        }
        returnFiles.push_back(peakPosFile);
 
        // 2. Write Interpolated Tangent as standard mode.dat
        AtomMatrix peakMode = (1.0 - f) * (*neb->tangent[leftIdx]) +
                              f * (*neb->tangent[leftIdx + 1]);
        peakMode.normalize();
 
        std::string peakModeFile =
            std::format("peak{:02d}_mode.dat", peakCount);
        {
          std::ofstream modeOut(peakModeFile);
          if (modeOut) {
            for (long row = 0; row < peakMode.rows(); ++row) {
              modeOut << std::format("{:.17g} {:.17g} {:.17g}\n",
                                     peakMode(row, 0), peakMode(row, 1),
                                     peakMode(row, 2));
            }
            returnFiles.push_back(peakModeFile);
          }
        }
 
        QUILL_LOG_INFO(
            m_log,
            "Generated MMF peak {:02d} at position {:.3f} (Energy: {:.3f} eV)",
            peakCount, posFraction, relativeEnergy);
        peakCount++;
      }
    }
  }
 
  returnFiles.push_back("neb.dat");
  neb->printImageData(true, std::numeric_limits<size_t>::max());
}
 
void NudgedElasticBandJob::printEndState(NudgedElasticBand::NEBStatus status) {
  QUILL_LOG_DEBUG(m_log, "Final state: ");
  if (status == NudgedElasticBand::NEBStatus::GOOD)
    QUILL_LOG_DEBUG(m_log, "Nudged elastic band, successful.");
  else if (status == NudgedElasticBand::NEBStatus::BAD_MAX_ITERATIONS)
    QUILL_LOG_DEBUG(m_log, "Nudged elastic band, too many iterations.");
  else
    QUILL_LOG_WARNING(m_log, "Unknown status: {}!", static_cast<int>(status));
  return;
}