NEBObjectiveFunction.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/NudgedElasticBand.h"
 
namespace eonc {
 
VectorXd NEBObjectiveFunction::getGradient(bool fdstep) {
  if (neb->movedAfterForceCall)
    neb->updateForces();
  const long seg = 3 * neb->atoms;
  VectorXd gradV(seg * neb->numImages);
  for (long i = 1; i <= neb->numImages; i++) {
    // Negate in-place during copy to avoid a second pass over 40KB
    gradV.segment(seg * (i - 1), seg) =
        -VectorXd::Map(neb->projectedForce[i]->data(), seg);
  }
  return gradV;
}
 
double NEBObjectiveFunction::getEnergy() {
  double Energy{0};
  for (long i = 1; i <= neb->numImages; i++) {
    Energy += neb->path[i]->getPotentialEnergy();
  }
  return Energy;
}
 
void NEBObjectiveFunction::setPositions(const VectorXd &x) {
  neb->movedAfterForceCall = true;
  for (long i = 1; i <= neb->numImages; i++) {
    neb->path[i]->setPositions(AtomMatrix::Map(
        x.segment(3 * neb->atoms * (i - 1), 3 * neb->atoms).data(), neb->atoms,
        3));
  }
}
 
VectorXd NEBObjectiveFunction::getPositions() {
  VectorXd posV;
  posV.resize(3 * neb->atoms * neb->numImages);
  for (long i = 1; i <= neb->numImages; i++) {
    posV.segment(3 * neb->atoms * (i - 1), 3 * neb->atoms) =
        VectorXd::Map(neb->path[i]->getPositions().data(), 3 * neb->atoms);
  }
  return posV;
}
 
int NEBObjectiveFunction::degreesOfFreedom() {
  return 3 * neb->numImages * neb->atoms;
}
 
bool NEBObjectiveFunction::isUncertain() {
  double maxMaxUnc = std::numeric_limits<double>::lowest();
  double currentMaxUnc{0};
  for (long idx = 0; idx <= neb->numImages; idx++) {
    currentMaxUnc = neb->path[idx]->getEnergyVariance();
    if (currentMaxUnc > maxMaxUnc) {
      maxMaxUnc = currentMaxUnc;
    }
  }
  bool unc_conv{maxMaxUnc > params.gp_surrogate_options.uncertainty};
  if (unc_conv) {
    this->status = NudgedElasticBand::NEBStatus::MAX_UNCERTAINTY;
  }
  return unc_conv;
}
 
bool NEBObjectiveFunction::isConverged() {
  bool force_conv = getConvergence() < params.neb_options.force_tolerance;
  return force_conv;
}
 
double NEBObjectiveFunction::getConvergence() {
  return neb->convergenceForce();
}
 
VectorXd NEBObjectiveFunction::difference(const VectorXd &a,
                                          const VectorXd &b) {
  VectorXd pbcDiff(3 * neb->numImages * neb->atoms);
  for (int i = 1; i <= neb->numImages; i++) {
    int n = (i - 1) * 3 * neb->atoms;
    int m = 3 * neb->atoms;
    pbcDiff.segment(n, m) =
        neb->path[i]->pbcV(a.segment(n, m) - b.segment(n, m));
  }
  return pbcDiff;
}
 
} // namespace eonc