DynamicsSaddleSearch.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/DynamicsSaddleSearch.h"
#include "eon/BondBoost.h"
#include "eon/Dynamics.h"
#include "eon/EigenmodeStrategy.h"
#include "eon/MinModeSaddleSearch.h"
#include "eon/NudgedElasticBand.h"
#include "eon/SafeMath.h"
 
#include <algorithm>
#include <cmath>
#include <filesystem>
#include <limits>
#include <stdexcept>
 
namespace eonc {
 
int DynamicsSaddleSearch::run() {
  std::vector<std::shared_ptr<Matter>> mdSnapshots;
  std::vector<double> mdTimes;
  QUILL_LOG_DEBUG(log, "Starting dynamics NEB saddle search");
 
  if (std::filesystem::exists("masses.dat")) {
    QUILL_LOG_DEBUG(log, "Found mass weights file");
    Eigen::VectorXd masses =
        eonc::helpers::loadMasses("masses.dat", saddle->numberOfAtoms());
    saddle->setMasses(masses);
    QUILL_LOG_DEBUG(log, "Applied mass weights");
  } else {
    QUILL_LOG_DEBUG(log, "No mass weights file found");
  }
 
  Dynamics dyn(saddle.get(), params);
  QUILL_LOG_DEBUG(
      log, "Initializing velocities from Maxwell-Boltzmann distribution");
  dyn.setTemperature(params.saddle_search_options().dynamics.temperature);
  dyn.setThermalVelocity();
 
  const double dt = params.dynamics_options().time_step;
  if (!(dt > 0.0)) {
    throw std::invalid_argument(
        "DynamicsSaddleSearch: time_step must be positive");
  }
  int dephaseSteps = static_cast<int>(
      std::floor(params.parallel_replica_options().dephase_time / dt + 0.5));
 
  while (true) {
 
    QUILL_LOG_DEBUG(log, "Dephasing: {} steps", dephaseSteps);
    // always start from the initial configuration
    *saddle = *reactant;
    dyn.setThermalVelocity();
 
    // Dephase MD trajectory
    for (int step = 1; step <= dephaseSteps; step++) {
      dyn.oneStep(step);
    }
 
    // Check to see if a transition occured
    Matter min(pot, params);
    min = *saddle;
    min.relax();
 
    if (min.compare(*reactant)) {
      QUILL_LOG_DEBUG(log, "Dephasing successful");
      break;
    } else {
      QUILL_LOG_DEBUG(log, "Transition occured during dephasing; Restarting");
      dephaseSteps /= 2;
      if (dephaseSteps < 1)
        dephaseSteps = 1;
    }
  }
 
  BondBoost bondBoost(saddle.get(), params);
  // setBiasPotential does not own this object. Clear it before bondBoost
  // leaves the stack, including on the early returns below.
  struct BiasGuard {
    Matter *matter{nullptr};
    ~BiasGuard() {
      if (matter == nullptr) {
        return;
      }
      matter->setBiasPotential(nullptr);
      matter->setBiasForces(AtomMatrix::Zero(matter->numberOfAtoms(), 3));
    }
  } biasGuard;
  if (params.hyperdynamics_options().bias_potential ==
      Hyperdynamics::BOND_BOOST) {
    QUILL_LOG_DEBUG(log, "Initializing Bond Boost");
    bondBoost.initialize();
    saddle->setBiasPotential(&bondBoost);
    biasGuard.matter = saddle.get();
  }
 
  int checkInterval = static_cast<int>(
      params.saddle_search_options().dynamics.state_check_interval /
          params.dynamics_options().time_step +
      0.5);
  // A zero or sub-step interval floors to 0. step % 0 is undefined, and a
  // state check shorter than one dynamics step still has to run.
  if (checkInterval < 1) {
    checkInterval = 1;
  }
  int recordInterval =
      static_cast<int>(params.saddle_search_options().dynamics.record_interval /
                           params.dynamics_options().time_step +
                       0.5);
 
  if (params.debug_options().write_movies) {
    if (!eonc::io::io_ok(saddle->matter2con("dynamics", false))) {
      QUILL_LOG_WARNING(log, "Failed to write dynamics movie header");
    }
  }
 
  for (int step = 1; step <= params.dynamics_options().steps; step++) {
    if (params.hyperdynamics_options().bias_potential ==
        Hyperdynamics::BOND_BOOST) {
      // oneStep() calls getAccelerations(), and therefore boost(), more
      // than once. Advance the rmd_time counter once per MD step.
      bondBoost.advance();
    }
    dyn.oneStep(step);
 
    if (recordInterval != 0 && step % recordInterval == 0) {
      QUILL_LOG_DEBUG(log, "recording configuration at step {} time {:.3f}",
                      step,
                      step * params.dynamics_options().time_step *
                          params.constants().timeUnit);
      // BUG FIX: was sharing ownership with saddle instead of copying
      auto snapshot = std::make_shared<Matter>(*saddle);
      mdSnapshots.push_back(snapshot);
      mdTimes.push_back(step * params.dynamics_options().time_step);
    }
 
    if (params.debug_options().write_movies) {
      if (!eonc::io::io_ok(saddle->matter2con("dynamics", true))) {
        QUILL_LOG_WARNING(log, "Failed to append dynamics movie frame");
      }
    }
 
    if (step % checkInterval == 0) {
      QUILL_LOG_DEBUG(log, "Minimizing trajectory, step {}", step);
 
      product = std::make_shared<Matter>(*saddle);
      product->relax(false, false);
 
      if (!product->compare(*reactant)) {
        QUILL_LOG_DEBUG(log, "Found new state");
        // A record interval that rounds to 0 stores nothing. refineTransition
        // then has no frame, and subscript 0 is outside the vector.
        int image = -1;
        if (!mdSnapshots.empty()) {
          image = refineTransition(mdSnapshots, product);
        }
        if (image < 0 || static_cast<size_t>(image) >= mdSnapshots.size() ||
            static_cast<size_t>(image) >= mdTimes.size()) {
          QUILL_LOG_DEBUG(log,
                          "No MD snapshots; using the detecting configuration");
          time = step * params.dynamics_options().time_step;
        } else {
          *saddle = *mdSnapshots[static_cast<size_t>(image)];
          QUILL_LOG_DEBUG(log, "Found transition at snapshot image {}", image);
          for (int ii = 0; ii < static_cast<int>(mdTimes.size()); ii++) {
            QUILL_LOG_DEBUG(log, "MDTimes[{}] = {:.3f}", ii,
                            mdTimes[ii] * params.constants().timeUnit);
          }
          // Subtract half the record interval to avoid systematic bias
          time = mdTimes[static_cast<size_t>(image)] -
                 params.saddle_search_options().dynamics.record_interval / 2.0;
          // Half a record interval before the leaving frame is the unbiased
          // crossing, and it must not fall before the preceding reactant frame.
          if (image > 0 && time < mdTimes[static_cast<size_t>(image - 1)]) {
            time = mdTimes[static_cast<size_t>(image - 1)];
          }
        }
        QUILL_LOG_DEBUG(log, "Transition time {:.2f} fs",
                        time * params.constants().timeUnit);
 
        NudgedElasticBand neb(reactant, product, params, pot);
 
        if (!params.saddle_search_options().dynamics.linear_interpolation) {
          QUILL_LOG_DEBUG(
              log, "Interpolating initial band through MD transition state");
          AtomMatrix reactantToSaddle =
              saddle->pbc(saddle->getPositions() - reactant->getPositions());
          AtomMatrix saddleToProduct =
              saddle->pbc(product->getPositions() - saddle->getPositions());
          QUILL_LOG_DEBUG(log, "Initial band saved to neb_initial_band.con");
          if (!eonc::io::io_ok(
                  neb.path[0]->matter2con("neb_initial_band.con", false))) {
            QUILL_LOG_WARNING(log, "Failed to write neb_initial_band.con");
          }
          int mid = neb.numImages / 2 + 1;
          for (int img = 1; img <= neb.numImages; img++) {
            if (img < mid) {
              double frac = static_cast<double>(img) / static_cast<double>(mid);
              neb.path[img]->setPositions(reactant->getPositions() +
                                          frac * reactantToSaddle);
            } else if (img > mid) {
              double frac = static_cast<double>(img - mid) /
                            static_cast<double>(neb.numImages - mid + 1);
              neb.path[img]->setPositions(saddle->getPositions() +
                                          frac * saddleToProduct);
            } else {
              neb.path[img]->setPositions(saddle->getPositions());
            }
            if (!eonc::io::io_ok(
                    neb.path[img]->matter2con("neb_initial_band.con", true))) {
              QUILL_LOG_WARNING(log, "Failed to append neb_initial_band frame");
            }
          }
          if (!eonc::io::io_ok(neb.path[neb.numImages + 1]->matter2con(
                  "neb_initial_band.con", true))) {
            QUILL_LOG_WARNING(log,
                              "Failed to append neb_initial_band endpoint");
          }
        } else {
          QUILL_LOG_DEBUG(
              log, "Linear interpolation between minima used for initial band");
          if (!eonc::io::io_ok(
                  neb.path[0]->matter2con("neb_initial_band.con", false))) {
            QUILL_LOG_WARNING(log, "Failed to write neb_initial_band.con");
          }
          for (int j = 1; j <= neb.numImages + 1; j++) {
            if (!eonc::io::io_ok(
                    neb.path[j]->matter2con("neb_initial_band.con", true))) {
              QUILL_LOG_WARNING(log, "Failed to append neb_initial_band frame");
            }
          }
        }
 
        AtomMatrix mode;
        if (params.neb_options().max_iterations > 0) {
          auto minModeMethod =
              eonc::buildEigenmodeStrategy(saddle, params, pot);
 
          neb.compute();
          neb.printImageData(true);
          int extremumImage = -1;
          int jExt = 0;
          for (jExt = 0; jExt < neb.numExtrema; jExt++) {
            if (neb.extremumCurvature[jExt] <
                params.saddle_search_options().dynamics.max_init_curvature) {
              extremumImage =
                  static_cast<int>(std::floor(neb.extremumPosition[jExt]));
              *saddle = *neb.path[extremumImage];
              double interpDist = neb.extremumPosition[jExt] -
                                  static_cast<double>(extremumImage);
              AtomMatrix bandDir =
                  saddle->pbc(neb.path[extremumImage + 1]->getPositions() -
                              neb.path[extremumImage]->getPositions());
              saddle->setPositions(interpDist * bandDir +
                                   saddle->getPositions());
              mode = saddle->pbc(neb.path[extremumImage + 1]->getPositions() -
                                 saddle->getPositions());
              eonc::safemath::safe_normalize_inplace(mode);
              eonc::eigenmodeCompute(*minModeMethod, saddle, mode);
              double ev = eonc::eigenmodeGetEigenvalue(*minModeMethod);
              QUILL_LOG_DEBUG(log, "extrema #{} has eigenvalue {:.8f}",
                              jExt + 1, ev);
 
              if (ev < 0) {
                QUILL_LOG_DEBUG(
                    log, "chose image {} (extrema #{}) as extremum image",
                    extremumImage, jExt + 1);
                break;
              } else {
                extremumImage = -1;
              }
            }
          }
 
          if (extremumImage != -1) {
            *saddle = *neb.path[extremumImage];
            double interpDist =
                neb.extremumPosition[jExt] - static_cast<double>(extremumImage);
            QUILL_LOG_DEBUG(log, "interpDistance {}", interpDist);
            AtomMatrix bandDir =
                saddle->pbc(neb.path[extremumImage + 1]->getPositions() -
                            neb.path[extremumImage]->getPositions());
            saddle->setPositions(interpDist * bandDir + saddle->getPositions());
            mode = saddle->pbc(neb.path[extremumImage + 1]->getPositions() -
                               saddle->getPositions());
            eonc::safemath::safe_normalize_inplace(mode);
          } else {
            QUILL_LOG_DEBUG(
                log, "no maxima found, using max energy non-endpoint image");
            double maxEnergy = -std::numeric_limits<double>::infinity();
            for (int img = 1; img <= neb.numImages; img++) {
              double U = neb.path[img]->getPotentialEnergy();
              if (U > maxEnergy) {
                maxEnergy = U;
                *saddle = *neb.path[img];
                mode = saddle->pbc(neb.path[img + 1]->getPositions() -
                                   saddle->getPositions());
                eonc::safemath::safe_normalize_inplace(mode);
              }
            }
            if (maxEnergy <= reactant->getPotentialEnergy()) {
              QUILL_LOG_DEBUG(log, "warning: no barrier found");
              return MinModeSaddleSearch::STATUS_BAD_NO_BARRIER;
            }
          }
        } else {
          // No NEB iterations: the middle image of the initial band is the
          // guess, and the dimer still needs an n-atom tangent. An empty
          // mode replaces the dimer direction and then indexes every atom.
          neb.maxEnergyImage = neb.numImages / 2 + 1;
          const int img = static_cast<int>(neb.maxEnergyImage);
          const int last = static_cast<int>(neb.path.size()) - 1;
          const int from = std::clamp(img, 0, last);
          const int to = std::clamp(img + 1, 0, last);
          mode = saddle->pbc(neb.path[to]->getPositions() -
                             neb.path[from]->getPositions());
          if (mode.norm() > 0.0) {
            mode.normalize();
          } else {
            mode = AtomMatrix::Zero(saddle->numberOfAtoms(), 3);
            if (mode.rows() > 0) {
              mode(0, 0) = 1.0;
            }
          }
        }
 
        QUILL_LOG_DEBUG(
            log, "Initial saddle guess saved to saddle_initial_guess.con");
        if (!eonc::io::io_ok(saddle->matter2con("saddle_initial_guess.con"))) {
          QUILL_LOG_WARNING(log, "Failed to write saddle_initial_guess.con");
        }
        MinModeSaddleSearch search = MinModeSaddleSearch(
            saddle, mode, reactant->getPotentialEnergy(), params, pot);
        int minModeStatus = search.run();
 
        if (minModeStatus != MinModeSaddleSearch::STATUS_GOOD) {
          QUILL_LOG_DEBUG(log, "error in min mode saddle search");
          return minModeStatus;
        }
 
        eigenvalue = search.getEigenvalue();
        eigenvector = search.getEigenvector();
        QUILL_LOG_DEBUG(log, "eigenvalue: {:.3f}", eigenvalue);
 
        double barrier =
            saddle->getPotentialEnergy() - reactant->getPotentialEnergy();
        QUILL_LOG_DEBUG(log, "found barrier of {:.3f}", barrier);
        mdSnapshots.clear();
        mdTimes.clear();
        return MinModeSaddleSearch::STATUS_GOOD;
      } else {
        QUILL_LOG_DEBUG(log, "Still in original state");
        mdTimes.clear();
        mdSnapshots.clear();
      }
    }
  }
 
  mdSnapshots.clear();
  time = params.dynamics_options().steps * params.dynamics_options().time_step;
  return MinModeSaddleSearch::STATUS_BAD_MD_TRAJECTORY_TOO_SHORT;
}

int DynamicsSaddleSearch::refineTransition(
    const std::vector<std::shared_ptr<Matter>> &snapshots,
    std::shared_ptr<Matter> prod) {
  int lo = 0;
  int hi = static_cast<int>(snapshots.size()) - 1;
  // One frame is index 0. Zero frames yield -1, which is not a subscript.
  if (hi <= 0) {
    return hi;
  }
 
  QUILL_LOG_DEBUG(log, "refining transition time");
 
  while ((hi - lo) > 1) {
    int mid = lo + (hi - lo) / 2;
    QUILL_LOG_DEBUG(log, "minimizing image {}", mid);
    Matter snap(pot, params);
    snap = *snapshots[mid];
    snap.relax(false);
 
    if (snap.compare(*reactant)) {
      QUILL_LOG_DEBUG(log, "image {} minimizes to reactant", mid);
      lo = mid;
    } else {
      QUILL_LOG_DEBUG(log, "image {} minimizes to product", mid);
      *prod = snap;
      hi = mid;
    }
  }
 
  // Adjacent brackets make (lo + hi) / 2 equal lo, the snapshot that
  // still minimizes to the reactant. The transition is the higher index.
  if (hi < 0) {
    return 0;
  }
  return hi;
}
 
double DynamicsSaddleSearch::getEigenvalue() { return eigenvalue; }
 
AtomMatrix DynamicsSaddleSearch::getEigenvector() { return eigenvector; }
 
} // namespace eonc