BondBoost.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/BondBoost.h"
#include "eon/EonLogger.h"
#include "eon/HelperFunctions.h"
 
#include <algorithm>
#include <cctype>
#include <cmath>
#include <ranges>
#include <stdexcept>
#include <string>
#include <unordered_set>
#include <vector>
 
namespace eonc {
 
const char Hyperdynamics::NONE[] = "none";
const char Hyperdynamics::BOND_BOOST[] = "bond_boost";
 
BondBoost::BondBoost(Matter *matt, const Parameters &params)
    : matter{matt},
      parameters{params} {
  if (!matter) {
    throw std::invalid_argument("BondBoost: null Matter");
  }
  nAtoms = matter->numberOfAtoms();
}
 
BondBoost::~BondBoost() = default;
 
void BondBoost::initialize() {
  nBBs = 0;
  nReg = 1;
 
  const std::string &balString =
      parameters.hyperdynamics_options().boost_atom_list;
  std::string lowered = balString;
  std::ranges::transform(lowered, lowered.begin(), [](unsigned char c) {
    return static_cast<char>(std::tolower(c));
  });
 
  BAList.clear();
  if (lowered.empty() || lowered == "all") {
    for (long i = 0; i < nAtoms; ++i) {
      if (!matter->getFixed(i)) {
        BAList.push_back(i);
      }
    }
  } else {
    const auto atoms = eonc::helpers::split_string_int(balString, ",");
    if (atoms.empty()) {
      throw std::invalid_argument(
          "hyperdynamics.boost_atom_list must be 'all' or a comma list of "
          "CON file-order indices");
    }
    std::unordered_set<long> seen;
    for (int raw : atoms) {
      const long row = matter->mapFileRow(static_cast<long>(raw));
      if (row < 0 || row >= nAtoms) {
        throw std::out_of_range(
            "hyperdynamics.boost_atom_list index out of range");
      }
      if (matter->getFixed(row)) {
        continue;
      }
      if (seen.insert(row).second) {
        BAList.push_back(row);
      }
    }
  }
  if (BAList.empty()) {
    throw std::runtime_error("BondBoost: no boostable atoms");
  }
 
  nBAs = static_cast<long>(BAList.size());
  const std::unordered_set<long> boosted(BAList.begin(), BAList.end());
  RAList.clear();
  RAList.reserve(static_cast<size_t>(nAtoms - nBAs));
  for (long i = 0; i < nAtoms; ++i) {
    if (boosted.find(i) == boosted.end()) {
      RAList.push_back(i);
    }
  }
  nRAs = static_cast<long>(RAList.size());
 
  nTABs = nBAs * (nBAs - 1) / 2 + nBAs * nRAs;
  TABAList.assign(static_cast<size_t>(2 * std::max(nTABs, 0L)), 0);
  TABLList.setZero(nTABs, 1);
  QUILL_LOG_DEBUG(log, "BondBoost: {} boost atoms, {} rest, {} tagged bonds",
                  nBAs, nRAs, nTABs);
}
 
long BondBoost::rmdSteps() const {
  const double dt = parameters.dynamics_options().time_step;
  if (!(dt > 0.0)) {
    return 0;
  }
  return static_cast<long>(parameters.hyperdynamics_options().rmd_time / dt);
}
 
void BondBoost::advance() {
  const long RMDS = rmdSteps();
  // nReg starts at 1, so a zero sample count never enters the average below
  // and never divides by RMDS. TABLList would stay the zeros from
  // initialize(), and BondSelect would keep every tagged pair at length 0.
  if (RMDS <= 0) {
    if (nBBs == 0) {
      TABLList = Rmdsteps();
      nBBs = BondSelect();
    }
    nReg++;
    return;
  }
  if (nReg <= RMDS) {
    // Equilibration: sample bond lengths once per MD step.
    Matrix<double, Eigen::Dynamic, 1> TABL_tmp = Rmdsteps();
    TABLList = TABLList + (1.0 / RMDS) * TABL_tmp;
    nReg++;
    if (nReg == RMDS + 1) {
      nBBs = BondSelect();
    }
    return;
  }
  if (nReg == RMDS + 1) {
    nBBs = BondSelect();
  }
  nReg++;
}
 
double BondBoost::boost() {
  const long RMDS = rmdSteps();
  if (RMDS > 0 && nReg <= RMDS) {
    return 0.0;
  }
  // Callers that never advance() still select on the first evaluation.
  // With no equilibration samples that selection has to measure lengths
  // first; SafeHyper and ParallelReplica advance() before they get here.
  if (nBBs == 0) {
    if (RMDS <= 0) {
      TABLList = Rmdsteps();
    }
    nBBs = BondSelect();
  }
  Epsr_Q.resize(nBBs);
  CBBLList.setZero(nBBs, 1);
  const double biasPot = Booststeps();
  Epsr_Q.clear();
  return biasPot;
}

double BondBoost::Booststeps() {
  const double QRR = parameters.hyperdynamics_options().qrr;
  const double PRR = parameters.hyperdynamics_options().prr;
  const double DVMAX = parameters.hyperdynamics_options().dvmax;
  if (nBBs <= 0 || !(QRR > 0.0)) {
    return 0.0;
  }
  const double nBBsD = static_cast<double>(nBBs);
 
  AtomMatrix addForces(nBBs, 3);
  AtomMatrix TADF(nAtoms, 3);
  addForces.setZero();
  TADF.setZero();
 
  // Measure current bond lengths
  for (long i = 0; i < nBBs; i++) {
    CBBLList(i, 0) = matter->distance(BBAList[2 * i], BBAList[2 * i + 1]);
  }
 
  // Compute strain parameters and find maximum
  double epsrMax = 0.0;
  for (long i = 0; i < nBBs; i++) {
    const double eq = EBBLList(i, 0);
    if (!(eq > 0.0)) {
      Epsr_Q[i] = 0.0;
      continue;
    }
    Epsr_Q[i] = (CBBLList(i, 0) - eq) / eq / QRR;
    if (std::abs(Epsr_Q[i]) >= epsrMax) {
      epsrMax = std::abs(Epsr_Q[i]);
    }
  }
 
  // Envelope function A(eps_max)
  double A_eps = (1.0 - epsrMax * epsrMax) * (1.0 - epsrMax * epsrMax) /
                 (1.0 - PRR * PRR * epsrMax * epsrMax);
 
  // Sum of individual bias potentials
  double sumV = 0.0;
  if (epsrMax < 1.0) {
    for (long i = 0; i < nBBs; i++) {
      sumV += DVMAX * (1.0 - Epsr_Q[i] * Epsr_Q[i]) / nBBsD;
    }
  } else {
    A_eps = 0.0;
  }
 
  double boostFact = A_eps * sumV;
 
  // Compute bias forces per bond, accumulate on atoms
  for (long i = 0; i < nBBs; i++) {
    double dforce = 0.0;
    double fact1 =
        2.0 * A_eps * DVMAX * Epsr_Q[i] / QRR / EBBLList(i, 0) / nBBsD;
 
    if (std::abs(Epsr_Q[i]) < epsrMax) {
      dforce = fact1;
    } else {
      // Bond at maximum strain: additional envelope derivative
      double fTmp1 = 1.0 - PRR * PRR * Epsr_Q[i] * Epsr_Q[i];
      double fTmp2 = 1.0 - Epsr_Q[i] * Epsr_Q[i];
      double fact2 = 2.0 * fTmp2 * Epsr_Q[i] *
                     (2.0 * fTmp1 - PRR * PRR * fTmp2) / QRR / EBBLList(i, 0) /
                     fTmp1 / fTmp1;
      dforce = fact1 + sumV * fact2;
    }
 
    long a1 = BBAList[2 * i];
    long a2 = BBAList[2 * i + 1];
    if (!(EBBLList(i, 0) > 0.0)) {
      continue;
    }
    // distance() minimum-images the whole bond. pdistance minimum-images one
    // Cartesian axis after zeroing the other two, which is not that vector
    // when the cell is non-orthogonal.
    AtomMatrix delta(1, 3);
    delta.row(0) =
        matter->getPositions().row(a1) - matter->getPositions().row(a2);
    delta = matter->pbc(delta);
    const double R = delta.norm();
    if (!(R > 0.0)) {
      continue;
    }
 
    for (int j = 0; j < 3; j++) {
      const double rij = delta(0, j);
      const double fij = rij / R * dforce;
      addForces(i, j) = fij;
      TADF(a1, j) += fij;
      TADF(a2, j) -= fij;
    }
  }
 
  // Apply free-atom mask and set bias forces
  AtomMatrix biasForces = TADF.array() * matter->getFree().array();
  matter->setBiasForces(biasForces);
  return boostFact;
}

Matrix<double, Eigen::Dynamic, 1> BondBoost::Rmdsteps() {
  Matrix<double, Eigen::Dynamic, 1> bondLengths(nTABs, 1);
  long count = 0;
 
  // Boost-atom pairs
  for (long i = 0; i < nBAs; i++) {
    for (long j = i + 1; j < nBAs; j++) {
      bondLengths(count, 0) = matter->distance(BAList[i], BAList[j]);
      TABAList[2 * count] = BAList[i];
      TABAList[2 * count + 1] = BAList[j];
      count++;
    }
  }
 
  // Boost-atom to rest-atom pairs
  for (long i = 0; i < nBAs; i++) {
    for (long j = 0; j < nRAs; j++) {
      bondLengths(count, 0) = matter->distance(BAList[i], RAList[j]);
      TABAList[2 * count] = BAList[i];
      TABAList[2 * count + 1] = RAList[j];
      count++;
    }
  }
 
  if (count != nTABs) {
    QUILL_LOG_DEBUG(log, "Total involved bond count does not match expected\n");
  }
  return bondLengths;
}

long BondBoost::BondSelect() {
  const double qCutoff = parameters.hyperdynamics_options().qcut;
 
  // Count bonds within cutoff
  long nSelected = 0;
  for (long i = 0; i < nTABs; i++) {
    if (TABLList(i, 0) <= qCutoff) {
      nSelected++;
    }
  }
 
  EBBLList.setZero(nSelected, 1);
  BBAList.resize(2 * nSelected);
  long count = 0;
  for (long i = 0; i < nTABs; i++) {
    if (TABLList(i, 0) <= qCutoff) {
      EBBLList(count, 0) = TABLList(i, 0);
      BBAList[2 * count] = TABAList[2 * i];
      BBAList[2 * count + 1] = TABAList[2 * i + 1];
      count++;
    }
  }
  return nSelected;
}
 
} // namespace eonc