ConjugateGradients.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/ConjugateGradients.h"
#include "eon/SafeMath.h"
 
#include <cmath>
 
Eigen::VectorXd ConjugateGradients::getStep() {
  double a = std::fabs(m_force.dot(m_forceOld));
  double b = m_forceOld.squaredNorm();
  double gamma = 0.0;
  if (a < 0.5 * b) {
    // Polak-Ribiere way to determine how much to mix in of old direction
    gamma = eonc::safemath::safe_div(m_force.dot(m_force - m_forceOld), b, 0.0);
  } else {
    gamma = 0;
  }
  m_direction = m_force + gamma * m_directionOld;
  m_directionNorm = m_direction;
  eonc::safemath::safe_normalize_inplace(m_directionNorm);
  m_directionOld = m_direction;
  m_forceOld = m_force;
 
  // Only if value for max nr of iteration before reset
  if (m_optConfig.opts.cg.max_iter_before_reset > 0 &&
      m_optConfig.opts.cg.max_iter_before_reset <= m_cg_i) {
    m_cg_i = 0;
    m_forceOld.setZero();
    m_directionOld.setZero();
  }
  m_cg_i += 1;
 
  return m_direction;
}
 
int ConjugateGradients::step(double a_maxMove) {
  bool converged;
  if (m_optConfig.opts.cg.line_search) {
    converged = line_search(a_maxMove);
  } else {
    converged = single_step(a_maxMove);
  }
  if (converged)
    return 1;
  return 0;
}
 
int ConjugateGradients::line_search(double a_maxMove) {
  Eigen::VectorXd pos;
  Eigen::VectorXd posStep;
  Eigen::VectorXd forceBeforeStep;
  double stepSize;
  double projectedForce;
  double projectedForceBeforeStep;
  double curvature;
 
  forceBeforeStep = -m_objf->getGradient();
  m_force = forceBeforeStep;
  getStep();
  pos = m_objf->getPositions();
  projectedForceBeforeStep = m_force.dot(m_directionNorm);
 
  // move system an infinitesimal step to determine the optimal step size along
  // the search line
  posStep = pos + m_directionNorm * m_optConfig.finiteDifference;
  m_objf->setPositions(posStep);
  m_force = -m_objf->getGradient(true);
  projectedForce = m_force.dot(m_directionNorm);
  stepSize = m_optConfig.finiteDifference;
 
  int line_i = 0;
  do {
    // Determine curvature from last step (Secant method)
    curvature = std::fabs(eonc::safemath::safe_div(
        projectedForceBeforeStep - projectedForce, stepSize, 0.0));
    stepSize = eonc::safemath::safe_div(projectedForce, curvature, a_maxMove);
 
    if (a_maxMove < std::fabs(stepSize)) {
      // first part get the sign of stepSize
      stepSize = ((stepSize > 0) - (stepSize < 0)) * a_maxMove;
    }
 
    forceBeforeStep = m_force;
    projectedForceBeforeStep = projectedForce;
 
    pos += stepSize * m_directionNorm;
    m_objf->setPositions(pos);
    m_force = -m_objf->getGradient();
    projectedForce = m_force.dot(m_directionNorm);
 
    line_i += 1;
 
    // Line search considered converged based in the ratio between the projected
    // force and the norm of the true force
  } while (m_optConfig.opts.cg.line_converged <
               std::fabs(projectedForce) /
                   (std::sqrt(m_force.dot(m_force) +
                              m_optConfig.opts.cg.line_converged)) &&
           line_i < m_optConfig.opts.cg.line_search_max_iter);
  if (m_objf->isConverged())
    return 1;
  return 0;
}
 
int ConjugateGradients::single_step(double a_maxMove) {
  Eigen::VectorXd pos;
  Eigen::VectorXd posStep;
  Eigen::VectorXd forceAfterStep;
 
  m_force = -m_objf->getGradient();
  pos = m_objf->getPositions();
  getStep();
 
  // move system an infinitesimal step to determine the optimal step size along
  // the search line
  posStep = pos + m_directionNorm * m_optConfig.finiteDifference;
  m_objf->setPositions(posStep);
  forceAfterStep = -m_objf->getGradient(true);
 
  // Determine curvature
  double projectedForce1 = m_force.dot(m_directionNorm);
  double projectedForce2 = forceAfterStep.dot(m_directionNorm);
  double curvature = eonc::safemath::safe_div(
      projectedForce1 - projectedForce2, m_optConfig.finiteDifference, 0.0);
 
  double stepSize = a_maxMove;
 
  if (curvature > 0.0) {
    stepSize = eonc::safemath::safe_div(projectedForce1, curvature, a_maxMove);
  }
 
  if (m_optConfig.bowlBreakout && a_maxMove < 0.0) {
    stepSize = -a_maxMove;
    a_maxMove = -a_maxMove;
  }
 
  if (!m_optConfig.opts.cg.no_overshooting) {
    if (m_optConfig.bowlBreakout) {
      // max displacement is based on system not single atom
      pos += eonc::helpers::maxMotionAppliedV(stepSize * m_directionNorm,
                                              a_maxMove);
    } else {
      pos += eonc::helpers::maxAtomMotionAppliedV(stepSize * m_directionNorm,
                                                  a_maxMove);
    }
    m_objf->setPositions(pos);
  } else {
    // negative if product of the projected forces before and after the step are
    // in opposite directions
    double passedMinimum = -1.;
    double forceChange = 0.;
    while (passedMinimum < 0.0 &&
           0.1 * std::fabs(projectedForce1) < std::fabs(projectedForce2)) {
      posStep = pos + eonc::helpers::maxAtomMotionAppliedV(
                          stepSize * m_directionNorm, a_maxMove);
      m_objf->setPositions(posStep);
      forceAfterStep = -m_objf->getGradient(true);
      projectedForce2 = forceAfterStep.dot(m_directionNorm);
 
      passedMinimum = projectedForce1 * projectedForce2;
      if (passedMinimum < 0.0 &&
          0.1 * std::fabs(projectedForce1) < std::fabs(projectedForce2)) {
        forceChange = (projectedForce1 - projectedForce2);
        stepSize = eonc::safemath::safe_div(projectedForce1, forceChange, 0.0) *
                   stepSize;
        QUILL_LOG_DEBUG(m_log, "Force changed {}, step size adjusted to {}",
                        forceChange, stepSize);
      }
    }
  }
  if (m_optConfig.opts.cg.knock_out_max_move) {
    if (stepSize >= a_maxMove) {
      // knockout old search direction
      m_directionOld.setZero();
      m_forceOld.setZero();
      QUILL_LOG_DEBUG(m_log, "Resetting the old search direction");
    }
  }
 
  return m_objf->isConverged() ? 1 : 0;
}
 
int ConjugateGradients::run(size_t a_maxIterations, double a_maxMove) {
  size_t iterations = 0;
  while (!m_objf->isConverged() && iterations <= a_maxIterations) {
    step(a_maxMove);
    iterations++;
  }
  return m_objf->isConverged() ? 1 : 0;
}