Prefactor.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/Prefactor.h"
#include "eon/HelperFunctions.h"
#include "eon/Hessian.h"
 
#include "eon/EonLogger.h"
#include <algorithm>
#include <cmath>
#include <format>
#include <fstream>
 
int eonc::Prefactor::getPrefactors(const Parameters &parameters, Matter *min1,
                                   Matter *saddle, Matter *min2, double &pref1,
                                   double &pref2) {
  VectorXd min1Freqs, saddleFreqs, min2Freqs;
 
  // determine which atoms moved in the process
  VectorXi atoms;
 
  if (parameters.prefactor_options.filter_scheme ==
      eonc::Prefactor::FILTER_FRACTION) {
    atoms = movedAtomsPct(parameters, min1, saddle, min2);
  } else {
    atoms = movedAtoms(parameters, min1, saddle, min2);
  }
 
  assert(atoms.rows() > 0);
 
  // calculate min1 frequencies
  Hessian hessian(parameters, min1);
  min1Freqs = hessian.getFreqs(min1, atoms);
  if (min1Freqs.size() == 0) {
    EONC_LOG_ERROR("[Prefactor] Bad hessian: min1");
    return -1;
  }
  // remove zero modes
  if (parameters.structure_comparison_options.check_rotation) {
    min1Freqs = hessian.removeZeroFreqs(min1Freqs);
  }
 
  // calculate saddle frequencies
  saddleFreqs = hessian.getFreqs(saddle, atoms);
  if (saddleFreqs.size() == 0) {
    EONC_LOG_ERROR("[Prefactor] Bad hessian: saddle");
    return -1;
  }
  // remove zero modes
  if (parameters.structure_comparison_options.check_rotation) {
    saddleFreqs = hessian.removeZeroFreqs(saddleFreqs);
  }
 
  // calculate min2 frequencies
  min2Freqs = hessian.getFreqs(min2, atoms);
  if (min2Freqs.size() == 0) {
    if (!parameters.main_options.quiet) {
      EONC_LOG_ERROR("[Prefactor] Bad hessian: min2");
    }
    return -1;
  }
  // remove zero modes
  if (parameters.structure_comparison_options.check_rotation) {
    min2Freqs = hessian.removeZeroFreqs(min2Freqs);
  }
 
  // check Hessian sizes
  // BUG FIX: original checked min1 vs saddle twice, never min2
  if ((min1Freqs.size() != saddleFreqs.size()) ||
      (min2Freqs.size() != saddleFreqs.size())) {
    if (!parameters.main_options.quiet) {
      EONC_LOG_ERROR("[Prefactor] Bad prefactor: Hessian sizes do not match");
    }
    return -1;
  }
 
  logFreqs(min1Freqs, "minimum 1");
  logFreqs(saddleFreqs, "saddle");
  logFreqs(min2Freqs, "minimum 2");
 
  // check for correct number of negative modes
  // Bound each loop by its own vector: removeZeroFreqs returns a shorter
  // vector than the 3N it was handed, and each of the three is shrunk by a
  // separate call.
  int numNegFreq = 0;
  for (int i = 0; i < min1Freqs.size(); i++) {
    if (min1Freqs(i) < 0) {
      EONC_LOG_DEBUG("[Prefactor] min1 had negative mode of {}", min1Freqs(i));
      numNegFreq++;
    }
  }
  if (numNegFreq != 0) {
    EONC_LOG_DEBUG("[Prefactor] Error: {} negative modes at min1", numNegFreq);
    return -1;
  }
 
  numNegFreq = 0;
  for (int i = 0; i < saddleFreqs.size(); i++) {
    if (saddleFreqs(i) < 0) {
      numNegFreq++;
    }
  }
  if (numNegFreq != 1) {
    EONC_LOG_DEBUG("Error: {} negative modes at saddle", numNegFreq);
    return -1;
  }
 
  numNegFreq = 0;
  for (int i = 0; i < min2Freqs.size(); i++) {
    if (min2Freqs(i) < 0) {
      numNegFreq++;
    }
  }
  if (numNegFreq != 0) {
    EONC_LOG_DEBUG("Error: {} negative modes at min2", numNegFreq);
    return -1;
  }
 
  // calculate the prefactors
  pref1 = 1.0;
  pref2 = 1.0;
 
  if (parameters.prefactor_options.rate == eonc::Prefactor::RATE_HTST) {
 
    // products are calculated this way in order to avoid overflow
    for (int i = 0; i < saddleFreqs.size(); i++) {
      pref1 *= min1Freqs[i];
      pref2 *= min2Freqs[i];
      if (saddleFreqs[i] > 0) {
        pref1 /= saddleFreqs[i];
        pref2 /= saddleFreqs[i];
      }
    }
    pref1 = std::sqrt(pref1) / (2 * eonc::helpers::pi * 10.18e-15);
    pref2 = std::sqrt(pref2) / (2 * eonc::helpers::pi * 10.18e-15);
  } else if (parameters.prefactor_options.rate ==
             eonc::Prefactor::RATE_QQHTST) {
    double kB_T = parameters.main_options.temperature * 8.617332e-5; // eV
    double h_bar = 6.582119e-16;                                     // eV*s
    double h = 4.135667e-15;                                         // eV*s
    double temp = (h_bar / (2.0 * kB_T));
 
    for (int i = 0; i < min1Freqs.size(); i++) {
      pref1 = pref1 * (std::sinh(temp * (std::sqrt(min1Freqs[i]) / 10.18e-15)));
      pref2 = pref2 * (std::sinh(temp * (std::sqrt(min2Freqs[i]) / 10.18e-15)));
 
      if (saddleFreqs[i] > 0) {
        pref1 =
            pref1 / (std::sinh(temp * (std::sqrt(saddleFreqs[i]) / 10.18e-15)));
        pref2 =
            pref2 / (std::sinh(temp * (std::sqrt(saddleFreqs[i]) / 10.18e-15)));
      }
    }
    pref1 = 2. * kB_T / (h)*pref1;
    pref2 = 2. * kB_T / (h)*pref2;
  }
  EONC_LOG_DEBUG("[Prefactor] reactant to product prefactor: {:.3e}", pref1);
  EONC_LOG_DEBUG("[Prefactor] product to reactant prefactor: {:.3e}", pref2);
  return 0;
}
 
void eonc::Prefactor::logFreqs(const VectorXd &freqs, const char *name) {
  std::ofstream outFile("freqs.dat", std::ios::app);
  if (!outFile.is_open())
    return;
 
  outFile << std::format("[Prefactor] Frequencies at {}\n", name);
  for (int i = 0; i < freqs.size(); i++) {
    outFile << "\n" << std::format("{:10.6f} ", freqs(i));
    if ((i + 1) % 5 == 0) {
      outFile << "\n";
    }
  }
  outFile << "\n";
}
 
VectorXi eonc::Prefactor::movedAtoms(const Parameters &parameters, Matter *min1,
                                     Matter *saddle, Matter *min2) {
  long nAtoms = saddle->numberOfAtoms();
 
  VectorXi moved(nAtoms);
  moved.setConstant(-1);
 
  AtomMatrix diffMin1 =
      saddle->pbc(saddle->getPositions() - min1->getPositions());
  AtomMatrix diffMin2 =
      saddle->pbc(saddle->getPositions() - min2->getPositions());
 
  diffMin1.array() *= saddle->getFree().array();
  diffMin2.array() *= saddle->getFree().array();
 
  int nMoved = 0;
  for (int i = 0; i < nAtoms; i++) {
    if ((diffMin1.row(i).norm() >
         parameters.prefactor_options.min_displacement) ||
        (diffMin2.row(i).norm() >
         parameters.prefactor_options.min_displacement)) {
      // Avoid Eigen Array == scalar in boolean context (C++20 / Eigen 3.3
      // manylinux toolchains reject Array comparison as bool).
      if (std::find(moved.data(), moved.data() + nMoved, i) ==
          moved.data() + nMoved) {
        moved[nMoved] = i;
        nMoved++;
      }
      for (int j = 0; j < nAtoms; j++) {
        double diffRSaddle = saddle->distance(i, j);
 
        if (diffRSaddle < parameters.prefactor_options.within_radius &&
            (!saddle->getFixed(j))) {
          if (std::find(moved.data(), moved.data() + nMoved, j) ==
              moved.data() + nMoved) {
            moved[nMoved] = j;
            nMoved++;
          }
        }
      }
    }
  }
  return moved.head(nMoved);
}
 
VectorXi eonc::Prefactor::movedAtomsPct(const Parameters &parameters,
                                        Matter *min1, Matter *saddle,
                                        Matter *min2) {
  long nAtoms = saddle->numberOfAtoms();
  long nFree = saddle->numberOfFreeAtoms();
 
  VectorXi moved(nAtoms);
  moved.setConstant(-1);
 
  AtomMatrix diffMin1 =
      saddle->pbc(saddle->getPositions() - min1->getPositions());
  AtomMatrix diffMin2 =
      saddle->pbc(saddle->getPositions() - min2->getPositions());
 
  diffMin1.array() *= saddle->getFree().array();
  diffMin2.array() *= saddle->getFree().array();
 
  VectorXd diff(nAtoms);
  diff.setConstant(0.0);
 
  QUILL_LOG_DEBUG(
      eonc::log::get(),
      "[Prefactor] including all atoms that make up {:.3f}% of the motion",
      parameters.prefactor_options.filter_fraction * 100);
  double sum = 0.0;
  int mini = 0;
  for (int i = 0; i < nAtoms; i++) {
    diff[i] = std::max(diffMin1.row(i).norm(), diffMin2.row(i).norm());
    sum += diff[i];
    if (diff[i] <= diff[mini]) {
      mini = i;
    }
  }
 
  EONC_LOG_DEBUG("[Prefactor] sum of atom distances moved {:.4f}", sum);
  EONC_LOG_DEBUG("[Prefactor] max moved atom distance: {:.4f}",
                 diff.maxCoeff());
 
  int nMoved = 0;
  double d = 0.0;
  while (d / sum <= parameters.prefactor_options.filter_fraction &&
         nMoved < nFree) {
    int maxi = mini;
    for (int i = 0; i < nAtoms; i++) {
      if (diff[i] >= diff[maxi]) {
        if (std::find(moved.data(), moved.data() + nMoved, i) ==
            moved.data() + nMoved) {
          maxi = i;
        }
      }
    }
    moved[nMoved] = maxi;
    nMoved++;
    d += diff[maxi];
  }
 
  int totalAtoms = nMoved;
  for (int i = 0; i < nMoved; i++) {
    for (int j = 0; j < nAtoms; j++) {
      if (moved[i] == j)
        continue;
 
      double diffRSaddle = saddle->distance(moved[i], j);
 
      if (diffRSaddle < parameters.prefactor_options.within_radius &&
          (!saddle->getFixed(j))) {
 
        if (std::find(moved.data(), moved.data() + totalAtoms, j) ==
            moved.data() + totalAtoms) {
          moved[totalAtoms++] = j;
        }
      }
    }
  }
  EONC_LOG_DEBUG("[Prefactor] including {} atoms in the hessian ({} moved + {} "
                 "neighbors)",
                 totalAtoms, nMoved, totalAtoms - nMoved);
  return moved.head(totalAtoms);
}
 
VectorXi eonc::Prefactor::allFreeAtoms(Matter *matter) {
  long nAtoms = matter->numberOfAtoms();
 
  VectorXi moved(nAtoms);
  moved.setConstant(-1);
 
  int nMoved = 0;
  for (int i = 0; i < nAtoms; i++) {
    if (!matter->getFixed(i)) {
      moved[nMoved] = i;
      nMoved++;
    }
  }
  return moved.head(nMoved);
}