FIRE.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/FIRE.h"
#include "eon/HelperFunctions.h"
 
#include <algorithm>
#include <format>
#include <stdexcept>
 
namespace eonc {
 
int FIRE::step(double a_maxMove) {
  double P = 0;
  // Check convergence.
  if (m_objf->isConverged()) {
    return 1;
  }
 
  // Velocity Verlet
  Eigen::VectorXd f = -m_objf->getGradient();
  Eigen::VectorXd x = m_objf->getPositions();
 
  m_vel += f * m_dt;
  Eigen::VectorXd dx = m_vel * m_dt;
 
  dx = eonc::geometry::maxAtomMotionAppliedV(dx, a_maxMove);
  m_objf->setPositions(x + dx);
 
  f = -m_objf->getGradient();
  const double fnorm = f.norm();
  Eigen::VectorXd f_unit = Eigen::VectorXd::Zero(f.size());
  if (fnorm > 0.0) {
    f_unit = f / fnorm;
  }
 
  // FIRE
  P = f.dot(m_vel);
  m_vel = (1 - m_alpha) * m_vel + m_alpha * f_unit * m_vel.norm();
  QUILL_LOG_DEBUG(m_log,
                  "P: {:.4f}, v: {:.4f}, m_dt: {:.4f}, m_alpha: {:.4f}, N: {}",
                  P, m_vel.norm(), m_dt, m_alpha, m_N);
  if (P >= 0) {
    m_N++;
    if (m_N > m_N_min) {
      m_dt = std::min(m_dt * m_f_inc, m_dt_max);
      m_alpha = m_alpha * m_f_a;
    }
  } else {
    m_dt = m_dt * m_f_dec;
    m_vel.setZero();
    m_alpha = m_alpha_start;
    m_N = 0;
  }
 
  // add a sanity check on m_dt
  if (m_dt < 1e-6) {
    QUILL_LOG_CRITICAL(m_log, "[FIRE] [critical] m_dt is too small: {:.4f}",
                       m_dt);
    throw std::runtime_error(
        std::format("[FIRE] m_dt is too small: {:.4f}", m_dt));
  }
 
  m_iteration++;
  return m_objf->isConverged() ? 1 : 0;
}
 
int FIRE::run(size_t a_maxIterations, double a_maxMove) {
  while (!m_objf->isConverged() && m_iteration < a_maxIterations) {
    step(a_maxMove);
  }
  return m_objf->isConverged() ? 1 : 0;
}
 
} // namespace eonc