NEBOcinebController.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/NEBOcinebController.h"
#include "eon/EonLogger.h"
#include "eon/MinModeSaddleSearch.h"
#include "eon/NudgedElasticBand.h"
#include "eon/eonExceptions.hpp"
#include <algorithm>
#include <cmath>
 
namespace eonc::neb {
 
OCINEBController::Config
OCINEBController::fromParams(const Parameters &params) {
  auto &ci = params.neb_options().climbing_image;
  auto &r = ci.ocineb;
  return Config{
      r.use_mmf,
      r.trigger_force,
      r.trigger_factor,
      r.max_steps,
      r.ci_stability_count,
      r.angle_tol,
      params.neb_options().force_tolerance,
      r.restore_unhelpful,
  };
}
 
OCINEBController::OCINEBController(const Config &cfg)
    : cfg_{cfg} {}
 
void OCINEBController::initBaseline(double baseline_force) {
  baseline_force_ = baseline_force;
  current_threshold_ = std::max(baseline_force_ * cfg_.trigger_factor,
                                2.0 * cfg_.force_tolerance);
}
 
bool OCINEBController::shouldTrigger(double convForce, bool ci_active,
                                     long climbingImage, long numImages,
                                     int ciStabilityCounter) const {
  if (!cfg_.use_mmf || !ci_active)
    return false;
  if (climbingImage <= 0 || climbingImage > numImages)
    return false;
  if (ciStabilityCounter <= static_cast<int>(cfg_.ci_stability_count))
    return false;
  if (convForce <= cfg_.force_tolerance)
    return false;
  return convForce < current_threshold_;
}
 
void OCINEBController::updateStability(long climbingImage) {
  if (climbingImage == previousClimbingImage_) {
    ciStabilityCounter_++;
  } else {
    ciStabilityCounter_ = 0;
    previousClimbingImage_ = climbingImage;
    has_cached_mode_ = false;
  }
}
 
OCINEBController::MMFResult OCINEBController::run(eonc::NudgedElasticBand &neb,
                                                  double convForce) {
  auto *log = eonc::log::get();
 
  QUILL_LOG_DEBUG(log,
                  "Triggering MMF.  Force: {:.4f}, Threshold: {:.4f} "
                  "({:.2f}x baseline)",
                  convForce, current_threshold_,
                  current_threshold_ / baseline_force_);
 
  // The dimer moves this image. updateForces() may then hop
  // climbingImage to a neighbor; band CI force is that neighbor.
  const long walked = neb.climbingImage;
  AtomMatrix savedPositions = neb.path[walked]->getPositions();
  // The saved geometry is evaluated already; a restore puts its forces
  // back instead of asking the potential again.
  const AtomMatrix savedForces = neb.path[walked]->getForcesRaw();
  const double savedEnergy = neb.path[walked]->getPotentialEnergy();
 
  double alignment = 0.0;
  int mmfResult = runDimer(neb, alignment);
 
  // Always update forces after MMF
  neb.movedAfterForceCall = true;
  neb.updateForces();
  double newForce = neb.convergenceForce();
 
  // A positive-curvature walk can put the band force under tolerance
  // while the climbing image sits in a well. Restore that image
  // before accepting the band.
  if (convergedClimb(newForce, cfg_.force_tolerance, mmfResult)) {
    QUILL_LOG_DEBUG(log, "NEB converged after MMF. Force: {:.4f}", newForce);
    return {newForce, true, false};
  }
 
  const double walkedForce = neb.path[walked]->getForcesFreeV().norm();
  const bool mmfHelped = walkHelped(walkedForce, convForce, mmfResult);
 
  if (mmfHelped) {
    updateThresholdSuccess(convForce, walkedForce);
    QUILL_LOG_DEBUG(log,
                    "MMF helped (status={}, walked={}, now_ci={}). "
                    "Walked force: {:.4f} -> {:.4f} ({:.2f}x baseline). "
                    "Band CI force: {:.4f}. New threshold: {:.4f}",
                    mmfResult, walked, neb.climbingImage, convForce,
                    walkedForce, walkedForce / baseline_force_, newForce,
                    current_threshold_);
  } else {
    // Frontiers OCI-NEB (doi:10.3389/fchem.2026.1807063, Algorithm 1)
    // restores only on positive curvature (status -2). Alignment
    // failure keeps the walked CI and applies the linear penalty.
    // restore_unhelpful also restores on alignment reject / force
    // increase on the walked image. Default false: the paper setting.
    const bool restore = cfg_.restore_unhelpful || mmfResult == -2;
    if (restore) {
      neb.path[walked]->setPositions(savedPositions);
      neb.path[walked]->setEvaluation(savedForces, savedEnergy);
      neb.movedAfterForceCall = true;
      has_cached_mode_ = false;
      newForce = convForce;
    }
    updateThresholdBackoff(alignment);
    QUILL_LOG_DEBUG(
        log,
        "MMF backoff (status={}, walked={}, now_ci={}). "
        "Walked force: {:.4f} -> {:.4f}, band CI: {:.4f}, "
        "Alignment:  {:.3f}. {}New threshold: {:.4f} ({:.2f}x baseline)",
        mmfResult, walked, neb.climbingImage, convForce, walkedForce, newForce,
        alignment, restore ? "Restored CI. " : "", current_threshold_,
        current_threshold_ / baseline_force_);
  }
 
  bool shouldReset =
      (savedPositions - neb.path[walked]->getPositions()).norm() >
      neb.params.optimizer_options().max_move *
          neb.params.neb_options().image_count;
 
  if (shouldReset) {
    QUILL_LOG_DEBUG(log, "Resetting optimization history.");
  }
 
  ciStabilityCounter_ = 0;
 
  return {newForce, false, shouldReset};
}
 
int OCINEBController::runDimer(eonc::NudgedElasticBand &neb,
                               double &alignment) {
  auto *log = eonc::log::get();
  alignment = 0.0;
 
  if (neb.climbingImage <= 0 || neb.climbingImage > neb.numImages) {
    QUILL_LOG_WARNING(log, "Invalid climbing image for MMF:  {}",
                      neb.climbingImage);
    return -1;
  }
 
  AtomMatrix initialMode;
  if (has_cached_mode_) {
    initialMode = cached_mode_;
  } else {
    initialMode = *neb.tangent[neb.climbingImage];
  }
  double tangentNorm = initialMode.norm();
  if (tangentNorm < 1e-8) {
    QUILL_LOG_WARNING(log, "Tangent too small for MMF initialization");
    return -1;
  }
  initialMode /= tangentNorm;
 
  auto tempMinModeSearch = std::make_shared<MinModeSaddleSearch>(
      neb.path[neb.climbingImage], initialMode,
      neb.path[neb.climbingImage]->getPotentialEnergy(), neb.params, neb.pot);
 
  int minModeStatus;
  try {
    minModeStatus = tempMinModeSearch->run(cfg_.max_steps);
  } catch (const eonc::DimerModeRestoredException &) {
    minModeStatus = MinModeSaddleSearch::STATUS_DIMER_RESTORED_BEST;
    QUILL_LOG_DEBUG(log, "MMF:  Dimer restored to best state");
  } catch (const eonc::DimerModeLostException &) {
    minModeStatus = MinModeSaddleSearch::STATUS_DIMER_LOST_MODE;
    QUILL_LOG_WARNING(log, "Dimer lost mode during MMF refinement");
  }
 
  mmf_iterations_used_ += tempMinModeSearch->iteration;
 
  double eigenvalue = tempMinModeSearch->getEigenvalue();
  if (eigenvalue > 0.0) {
    QUILL_LOG_WARNING(log,
                      "MMF skipped: Positive curvature detected (eig={:.4f}).",
                      eigenvalue);
    return -2;
  }
 
  AtomMatrix finalModeMatrix = tempMinModeSearch->getEigenvector();
  VectorXd finalMode =
      VectorXd::Map(finalModeMatrix.data(), finalModeMatrix.size());
  VectorXd currentTangent =
      VectorXd::Map(neb.tangent[neb.climbingImage]->data(),
                    neb.tangent[neb.climbingImage]->size());
  if (finalMode.size() != currentTangent.size()) {
    QUILL_LOG_WARNING(log,
                      "MMF mode size {} != tangent size {}; skip alignment",
                      finalMode.size(), currentTangent.size());
    alignment = 0.0;
  } else if (finalMode.norm() == 0.0 || currentTangent.norm() == 0.0) {
    alignment = 0.0;
  } else {
    alignment =
        std::abs(finalMode.normalized().dot(currentTangent.normalized()));
  }
 
  if (minModeStatus == MinModeSaddleSearch::STATUS_GOOD ||
      minModeStatus == MinModeSaddleSearch::STATUS_DIMER_RESTORED_BEST) {
    if (alignment < cfg_.angle_tol) {
      QUILL_LOG_WARNING(
          log,
          "MMF converged/restored but mode drifted (alignment={:.3f} < {:.3f})",
          alignment, cfg_.angle_tol);
      return -1;
    }
    cached_mode_ = finalModeMatrix;
    has_cached_mode_ = true;
    return 0;
  } else if (minModeStatus == MinModeSaddleSearch::STATUS_BAD_MAX_ITERATIONS) {
    return 1;
  } else {
    QUILL_LOG_WARNING(log, "MMF failed.  Mode-tangent alignment: {:.3f}",
                      alignment);
    return -1;
  }
}
 
void OCINEBController::updateThresholdSuccess(double convForce,
                                              double newForce) {
  current_threshold_ = newForce * (0.5 + 0.4 * (newForce / convForce));
  double max_threshold = std::max(baseline_force_ * cfg_.trigger_factor,
                                  2.0 * cfg_.force_tolerance);
  current_threshold_ = std::min(current_threshold_, max_threshold);
}
 
void OCINEBController::updateThresholdBackoff(double alignment) {
  double alpha = std::clamp(alignment, 0.0, 1.0);
  double penalty_factor = 0.5 + 0.5 * alpha;
  current_threshold_ = baseline_force_ * cfg_.trigger_factor * penalty_factor;
  // Lower bound on the MMF trigger threshold. Scaled by force_tolerance so
  // the MMF gate never collapses to zero when the NEB is already near
  // convergence, but also capped by the trigger_factor envelope so a
  // loose force_tolerance cannot push min_threshold above the cap and
  // starve MMF activation.
  double min_threshold = 2.0 * cfg_.force_tolerance;
  current_threshold_ = std::max(current_threshold_, min_threshold);
}
 
} // namespace eonc::neb