eonc::Hessian Class Reference

#include <Hessian.h>

Public Member Functions
	Hessian (const Parameters &params, Matter *matter)
	~Hessian ()=default
MatrixXd	getHessian (Matter *matterIn, const VectorXi &atomsIn)
VectorXd	getFreqs (Matter *matterIn, const VectorXi &atomsIn)
VectorXd	removeZeroFreqs (const VectorXd &freqs)
const MatrixXd &	getModes () const noexcept
	Eigenvectors of the mass-weighted Hessian, one column per eigenvalue of getFreqs(), over the mobile degrees of freedom.
void	writeHessianFile (bool on) noexcept
	Whether a finished Hessian goes to hessian.dat (on by default).

Private Member Functions
bool	calculate ()
bool	finalizeHessian (int size)
bool	calculateColored (double cutoff, double dr, FdScheme scheme)
bool	calculateSerial (double dr, FdScheme scheme)
bool	calculateBatched (double dr, FdScheme scheme)

Private Attributes
Matter *	matter
const Parameters &	parameters
MatrixXd	hessian
VectorXd	freqs
MatrixXd	modes
bool	writeFile = true
VectorXi	atoms
eonc::log::Scoped	log

Detailed Description

Definition at line 47 of file Hessian.h.

Constructor & Destructor Documentation

Hessian()

eonc::Hessian::Hessian	(	const Parameters &	params,
		Matter *	matter )

Definition at line 81 of file Hessian.cpp.

    : matter{matter},
      parameters{params} {
  hessian.resize(0, 0);
  freqs.resize(0);
}

~Hessian()

eonc::Hessian::~Hessian ( )

default

Member Function Documentation

calculate()

bool eonc::Hessian::calculate ( )

private

Definition at line 271 of file Hessian.cpp.

                        {
  int nAtoms = matter->numberOfAtoms();
 
  int size = static_cast<int>(atoms.rows()) * 3;
  QUILL_LOG_DEBUG(log, "[Hessian] Hessian size: {}\n", size);
  if (size == 0) {
    return false;
  }
 
  // Mobile-atom polarity: indices in `atoms` are FD-displaced DOF owners.
  for (int a = 0; a < atoms.rows(); ++a) {
    const long idx = atoms(a);
    if (idx < 0 || idx >= nAtoms) {
      QUILL_LOG_ERROR(log,
                      "[Hessian] atom index {} out of range [0, {}) at list "
                      "entry {}; aborting FD Hessian",
                      idx, nAtoms, a);
      return false;
    }
  }
 
  double dr = parameters.main_options().finiteDifference;
  if (!(dr > 0.0) || !std::isfinite(dr)) {
    QUILL_LOG_ERROR(log, "[Hessian] invalid finiteDifference dr={}\n", dr);
    return false;
  }
 
  const FdScheme scheme = parseFdScheme(parameters.hessian_options().fd_scheme);
  const std::string &ckptPath = parameters.hessian_options().checkpoint_path;
 
  hessian.resize(size, size);
  hessian.setZero();
 
  // Net-force removal adds the same shift to every atom, so columns are
  // no longer confined to the cutoff neighborhood. A column checkpoint is
  // also stored one coordinate at a time.
  bool anyFixed = false;
  for (int i = 0; i < nAtoms; ++i) {
    if (matter->getFixed(i)) {
      anyFixed = true;
      break;
    }
  }
  const bool netCoupled =
      parameters.main_options().removeNetForce && nAtoms > 1 && !anyFixed;
  double cutoff = 0.0;
  if (matter->getPotential()) {
    cutoff = matter->getPotential()->finiteCutoff();
  }
  if (!netCoupled && ckptPath.empty() && cutoff > 0.0 &&
      std::isfinite(cutoff)) {
    if (calculateColored(cutoff, dr, scheme)) {
      return true;
    }
    hessian.setZero();
  }
  // A potential that evaluates batches (calculator groups, GPU models)
  // takes the displaced structures together; the column checkpoint stays
  // on the one-column-at-a-time path.
  if (ckptPath.empty() && matter->getPotential() &&
      matter->getPotential()->supportsBatchEvaluation() && size > 1) {
    return calculateBatched(dr, scheme);
  }
  return calculateSerial(dr, scheme);
}

calculateBatched()

bool eonc::Hessian::calculateBatched ( double dr, FdScheme scheme )

private

Definition at line 337 of file Hessian.cpp.

                                                         {
  const int size = static_cast<int>(atoms.rows()) * 3;
  const AtomMatrix pos = matter->getPositions();
  auto pot = matter->getPotential();
 
  Matter base(*matter);
  const AtomMatrix force0 = base.getForces();
  if (!force0.allFinite()) {
    QUILL_LOG_ERROR(log, "[Hessian] non-finite forces at undisplaced geometry; "
                         "aborting FD Hessian");
    return false;
  }
 
  // Stencil points per column, in the order fdForceDerivative takes them.
  std::vector<double> steps{1.0};
  if (scheme != FdScheme::OneSided) {
    steps.push_back(-1.0);
  }
  if (scheme == FdScheme::Fourth) {
    steps.push_back(2.0);
    steps.push_back(-2.0);
  }
  const int perColumn = static_cast<int>(steps.size());
  const long nAtoms = matter->numberOfAtoms();
  const VectorXi nrs = matter->getAtomicNrs();
  const Matrix3d box =
      matter->getPeriodic() ? matter->getCell() : Matrix3d::Zero().eval();
 
  // Columns in chunks, so memory stays bounded for large mobile sets.
  constexpr int kChunkColumns = 32;
  std::vector<Matter> displaced;
  for (int c0 = 0; c0 < size; c0 += kChunkColumns) {
    const int c1 = std::min(size, c0 + kChunkColumns);
    const long n = static_cast<long>(c1 - c0) * perColumn;
    displaced.assign(static_cast<size_t>(n), base);
    std::vector<const double *> posPtr, boxPtr;
    std::vector<const int *> nrsPtr;
    std::vector<double *> frcPtr;
    for (int i = c0; i < c1; ++i) {
      for (int k = 0; k < perColumn; ++k) {
        Matter &m = displaced[static_cast<size_t>((i - c0) * perColumn + k)];
        AtomMatrix p = pos;
        p(atoms(i / 3), i % 3) += steps[static_cast<size_t>(k)] * dr;
        m.setPositions(p);
      }
    }
    for (auto &m : displaced) {
      posPtr.push_back(m.getPositions().data());
      nrsPtr.push_back(nrs.data());
      frcPtr.push_back(m.forcesData());
      boxPtr.push_back(box.data());
    }
    std::vector<double> energies(static_cast<size_t>(n)),
        variances(static_cast<size_t>(n));
    pot->forceBatch(n, nAtoms, posPtr.data(), nrsPtr.data(), frcPtr.data(),
                    energies.data(), variances.data(), boxPtr.data());
    for (long j = 0; j < n; ++j) {
      displaced[static_cast<size_t>(j)].setComputedPotential(
          energies[static_cast<size_t>(j)], variances[static_cast<size_t>(j)]);
    }
    for (int i = c0; i < c1; ++i) {
      auto forces = [&](int k) -> const AtomMatrix & {
        return displaced[static_cast<size_t>((i - c0) * perColumn + k)]
            .getForces();
      };
      const AtomMatrix &fPlus = forces(0);
      const AtomMatrix &fMinus = perColumn > 1 ? forces(1) : force0;
      const AtomMatrix &fPlus2 = perColumn > 2 ? forces(2) : force0;
      const AtomMatrix &fMinus2 = perColumn > 3 ? forces(3) : force0;
      if (!fPlus.allFinite() || !fMinus.allFinite() || !fPlus2.allFinite() ||
          !fMinus2.allFinite()) {
        QUILL_LOG_ERROR(log,
                        "[Hessian] non-finite forces for FD column {}; "
                        "aborting FD Hessian",
                        i);
        return false;
      }
      const AtomMatrix slope =
          fdForceDerivative(scheme, dr, force0, fPlus, fMinus, fPlus2, fMinus2);
      for (int j = 0; j < size; j++) {
        const double effMass = std::sqrt(matter->getMass(atoms(j / 3)) *
                                         matter->getMass(atoms(i / 3)));
        hessian(i, j) =
            eonc::safemath::safe_div(-slope(atoms(j / 3), j % 3), effMass, 0.0);
      }
    }
  }
  return finalizeHessian(size);
}

calculateColored()

bool eonc::Hessian::calculateColored ( double cutoff, double dr, FdScheme scheme )

private

Definition at line 427 of file Hessian.cpp.

                                                                        {
  const int nAtoms = static_cast<int>(matter->numberOfAtoms());
  const int nMobile = static_cast<int>(atoms.rows());
  const int size = nMobile * 3;
  const MobileColoring coloring = buildMobileColoring(*matter, atoms, cutoff);
  if (static_cast<int>(coloring.color.size()) != nMobile ||
      static_cast<int>(coloring.closed.size()) != nMobile) {
    return false;
  }
 
  int nColors = 0;
  for (int c : coloring.color) {
    if (c < 0) {
      return false;
    }
    nColors = std::max(nColors, c + 1);
  }
  if (nColors <= 0) {
    return false;
  }
  QUILL_LOG_DEBUG(log,
                  "[Hessian] cutoff coloring: {} colors for {} mobile atoms\n",
                  nColors, nMobile);
 
  std::vector<std::vector<int>> members(static_cast<std::size_t>(nColors));
  for (int ia = 0; ia < nMobile; ++ia) {
    members[static_cast<std::size_t>(
                coloring.color[static_cast<std::size_t>(ia)])]
        .push_back(ia);
  }
  std::vector<std::vector<int>> owner(
      static_cast<std::size_t>(nColors),
      std::vector<int>(static_cast<std::size_t>(nMobile), -1));
  for (int c = 0; c < nColors; ++c) {
    for (int ia : members[static_cast<std::size_t>(c)]) {
      for (int k : coloring.closed[static_cast<std::size_t>(ia)]) {
        int &slot =
            owner[static_cast<std::size_t>(c)][static_cast<std::size_t>(k)];
        if (slot >= 0 && slot != ia) {
          return false;
        }
        slot = ia;
      }
    }
  }
 
  Matter matterTemp(*matter);
  const AtomMatrix pos = matter->getPositions();
  AtomMatrix posDisplace(nAtoms, 3);
  AtomMatrix force0 = matterTemp.getForces();
  if (!force0.allFinite()) {
    QUILL_LOG_ERROR(log, "[Hessian] non-finite forces at undisplaced geometry; "
                         "aborting FD Hessian");
    return false;
  }
 
  for (int dir = 0; dir < 3; ++dir) {
    for (int c = 0; c < nColors; ++c) {
      const auto &group = members[static_cast<std::size_t>(c)];
      if (group.empty()) {
        continue;
      }
      auto shifted = [&](double scale, const char *side) -> AtomMatrix {
        posDisplace.setZero();
        for (int ia : group) {
          posDisplace(atoms(ia), dir) = scale * dr;
        }
        matterTemp.setPositions(pos + posDisplace);
        AtomMatrix force = matterTemp.getForces();
        if (!force.allFinite()) {
          QUILL_LOG_ERROR(log,
                          "[Hessian] non-finite forces for color {} dir {} "
                          "({}); aborting FD Hessian",
                          c, dir, side);
          return AtomMatrix();
        }
        return force;
      };
      const AtomMatrix forcePlus = shifted(1.0, "+");
      if (forcePlus.size() == 0) {
        return false;
      }
      AtomMatrix forceMinus = force0;
      AtomMatrix forcePlus2 = force0;
      AtomMatrix forceMinus2 = force0;
      if (scheme != FdScheme::OneSided) {
        forceMinus = shifted(-1.0, "-");
        if (forceMinus.size() == 0) {
          return false;
        }
      }
      if (scheme == FdScheme::Fourth) {
        forcePlus2 = shifted(2.0, "+2");
        if (forcePlus2.size() == 0) {
          return false;
        }
        forceMinus2 = shifted(-2.0, "-2");
        if (forceMinus2.size() == 0) {
          return false;
        }
      }
      const AtomMatrix slope = fdForceDerivative(
          scheme, dr, force0, forcePlus, forceMinus, forcePlus2, forceMinus2);
      const AtomMatrix zero = AtomMatrix::Zero(nAtoms, 3);
      for (int ia : group) {
        const int col = ia * 3 + dir;
        writeMassWeighted(hessian, *matter, atoms, col,
                          static_cast<int>(atoms(ia)), slope, zero, 1.0,
                          owner[static_cast<std::size_t>(c)], ia);
      }
    }
  }
  return finalizeHessian(size);
}

calculateSerial()

bool eonc::Hessian::calculateSerial ( double dr, FdScheme scheme )

private

Definition at line 542 of file Hessian.cpp.

                                                        {
  const int nAtoms = static_cast<int>(matter->numberOfAtoms());
  const int size = static_cast<int>(atoms.rows()) * 3;
  const std::string &ckptPath = parameters.hessian_options().checkpoint_path;
  const bool wantResume =
      parameters.hessian_options().resume && !ckptPath.empty();
 
  AtomMatrix pos = matter->getPositions();
  AtomMatrix posDisplace(nAtoms, 3);
  AtomMatrix posTemp(nAtoms, 3);
  AtomMatrix forcePlus(nAtoms, 3);
  AtomMatrix forceMinus(nAtoms, 3);
 
  int startCol = 0;
  if (wantResume && loadColumnCheckpoint(ckptPath, size, startCol, hessian)) {
    QUILL_LOG_DEBUG(log, "[Hessian] resume from column {} / {}\n", startCol,
                    size);
  } else {
    startCol = 0;
    hessian.setZero();
  }
 
  Matter matterTemp(*matter);
  AtomMatrix force0 = matterTemp.getForces();
  if (!force0.allFinite()) {
    QUILL_LOG_ERROR(log, "[Hessian] non-finite forces at undisplaced geometry; "
                         "aborting FD Hessian");
    return false;
  }
 
  for (int i = startCol; i < size; i++) {
    posDisplace.setZero();
    posDisplace(atoms(i / 3), i % 3) = dr;
 
    posTemp = pos + posDisplace;
    matterTemp.setPositions(posTemp);
    forcePlus = matterTemp.getForces();
    if (!forcePlus.allFinite()) {
      QUILL_LOG_ERROR(log,
                      "[Hessian] non-finite forces for FD column {} (+); "
                      "aborting FD Hessian",
                      i);
      return false;
    }
 
    if (scheme != FdScheme::OneSided) {
      posTemp = pos - posDisplace;
      matterTemp.setPositions(posTemp);
      forceMinus = matterTemp.getForces();
      if (!forceMinus.allFinite()) {
        QUILL_LOG_ERROR(log,
                        "[Hessian] non-finite forces for FD column {} (-); "
                        "aborting FD Hessian",
                        i);
        return false;
      }
    }
    AtomMatrix forcePlus2 = force0;
    AtomMatrix forceMinus2 = force0;
    if (scheme == FdScheme::Fourth) {
      posDisplace(atoms(i / 3), i % 3) = 2.0 * dr;
      matterTemp.setPositions(pos + posDisplace);
      forcePlus2 = matterTemp.getForces();
      if (!forcePlus2.allFinite()) {
        QUILL_LOG_ERROR(log,
                        "[Hessian] non-finite forces for FD column {} (+2); "
                        "aborting FD Hessian",
                        i);
        return false;
      }
      matterTemp.setPositions(pos - posDisplace);
      forceMinus2 = matterTemp.getForces();
      if (!forceMinus2.allFinite()) {
        QUILL_LOG_ERROR(log,
                        "[Hessian] non-finite forces for FD column {} (-2); "
                        "aborting FD Hessian",
                        i);
        return false;
      }
    }
    // H ≈ -dF, mass-weighted. dF is the selected real stencil.
    const AtomMatrix slope = fdForceDerivative(
        scheme, dr, force0, forcePlus, forceMinus, forcePlus2, forceMinus2);
    for (int j = 0; j < size; j++) {
      hessian(i, j) = -slope(atoms(j / 3), j % 3);
      const double effMass = std::sqrt(matter->getMass(atoms(j / 3)) *
                                       matter->getMass(atoms(i / 3)));
      hessian(i, j) = eonc::safemath::safe_div(hessian(i, j), effMass, 0.0);
    }
 
    if (!ckptPath.empty()) {
      // next column to compute after a clean interrupt
      saveColumnCheckpoint(ckptPath, size, i + 1, hessian);
    }
  }
  return finalizeHessian(size);
}

finalizeHessian()

bool eonc::Hessian::finalizeHessian ( int size )

private

Definition at line 640 of file Hessian.cpp.

                                      {
  const std::string &ckptPath = parameters.hessian_options().checkpoint_path;
 
  // Symmetrize (FD noise breaks H=H^T; required for vib analysis)
  for (int i = 0; i < size; i++) {
    for (int j = 0; j < i; j++) {
      hessian(i, j) = (hessian(i, j) + hessian(j, i)) / 2;
      hessian(j, i) = hessian(i, j);
    }
  }
 
  if (!hessian.allFinite()) {
    QUILL_LOG_ERROR(log, "[Hessian] non-finite entries after FD assembly; "
                         "aborting eigen solve");
    return false;
  }
 
  if (!parameters.main_options().quiet) {
    QUILL_LOG_DEBUG(log, "[Hessian] writing hessian\n");
  }
  if (writeFile) {
    // A previous case in the same process can still hold hessian.dat on
    // Windows. Replace it, and try once more after removing the old file.
    std::ofstream hessfile("hessian.dat", std::ios::out | std::ios::trunc);
    if (!hessfile) {
      std::remove("hessian.dat");
      hessfile.clear();
      hessfile.open("hessian.dat", std::ios::out | std::ios::trunc);
    }
    if (!hessfile) {
      QUILL_LOG_ERROR(log, "[Hessian] failed to open hessian.dat");
      return false;
    }
    hessfile << hessian;
    hessfile.close();
    if (!hessfile) {
      QUILL_LOG_ERROR(log, "[Hessian] failed to write hessian.dat");
      return false;
    }
  }
 
  // Completed run: remove checkpoint so a later job does not resume stale cols
  if (!ckptPath.empty()) {
    std::remove(ckptPath.c_str());
  }
 
  double t0, t1;
  eonc::helpers::getTime(&t0, nullptr, nullptr);
  QUILL_LOG_DEBUG(log, "[Hessian] calculating eigen values of the hessian\n");
  // ColMajor copy for SelfAdjointEigenSolver (eOn MatrixXd is RowMajor)
  using ColMajorXd =
      Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>;
  ColMajorXd hessianCol = hessian;
  const bool withModes = parameters.hessian_options().write_modes;
  Eigen::SelfAdjointEigenSolver<ColMajorXd> es(
      hessianCol,
      withModes ? Eigen::ComputeEigenvectors : Eigen::EigenvaluesOnly);
  eonc::helpers::getTime(&t1, nullptr, nullptr);
  QUILL_LOG_DEBUG(log, "[Hessian] eigenvalue problem took {:.4e} seconds\n",
                  t1 - t0);
  if (es.info() != Eigen::Success) {
    QUILL_LOG_ERROR(log,
                    "[Hessian] SelfAdjointEigenSolver failed (info={}); "
                    "aborting",
                    static_cast<int>(es.info()));
    return false;
  }
  freqs = es.eigenvalues();
  if (!freqs.allFinite()) {
    QUILL_LOG_ERROR(log, "[Hessian] non-finite eigenvalues; aborting");
    return false;
  }
  if (withModes) {
    modes = es.eigenvectors();
  } else {
    modes.resize(0, 0);
  }
 
  return true;
}

getFreqs()

VectorXd eonc::Hessian::getFreqs	(	Matter *	matterIn,
		const VectorXi &	atomsIn )

Definition at line 102 of file Hessian.cpp.

                                                                    {
  if ((matter != matterIn) || (atoms.size() != atomsIn.size()) ||
      (atoms != atomsIn) || (hessian.rows() == 0)) {
    hessian.resize(0, 0);
    matter = matterIn;
    atoms = atomsIn;
 
    if (!calculate()) {
      freqs.resize(0);
      hessian.resize(0, 0);
    }
  }
  return freqs;
}

getHessian()

MatrixXd eonc::Hessian::getHessian	(	Matter *	matterIn,
		const VectorXi &	atomsIn )

Definition at line 88 of file Hessian.cpp.

                                                                      {
  if ((matter != matterIn) || (atoms.size() != atomsIn.size()) ||
      (atoms != atomsIn) || (hessian.rows() == 0)) {
    hessian.resize(0, 0);
    matter = matterIn;
    atoms = atomsIn;
 
    if (!calculate()) {
      hessian.resize(0, 0);
    }
  }
  return hessian;
}

getModes()

const MatrixXd & eonc::Hessian::getModes ( ) const

inlinenodiscardnoexcept

Eigenvectors of the mass-weighted Hessian, one column per eigenvalue of getFreqs(), over the mobile degrees of freedom.

Empty unless [Hessian] write_modes is set.

Definition at line 58 of file Hessian.h.

{ return modes; }

removeZeroFreqs()

VectorXd eonc::Hessian::removeZeroFreqs

(

const VectorXd &

freqs

)

Definition at line 721 of file Hessian.cpp.

                                                       {
  QUILL_LOG_DEBUG(log, "[Hessian] removing zero frequency modes");
  int size = freqs.size();
  if (size != 3 * matter->numberOfAtoms()) {
    return freqs;
  }
  VectorXd newfreqs;
  newfreqs.resize(size);
  int nremoved = 0;
  for (int i = 0; i < size; i++) {
    if (std::abs(freqs(i)) > parameters.hessian_options().zero_freq_value) {
      newfreqs(i - nremoved) = freqs(i);
    } else {
      nremoved++;
    }
  }
 
  if (!trivialModeCountIsPhysical(nremoved, matter->numberOfFixedAtoms())) {
    QUILL_LOG_WARNING(log,
                      "[Hessian] found {} trivial eigenmodes; a free cluster "
                      "has 6 (5 if linear), a periodic cell 3, and a "
                      "structure with fixed atoms none",
                      nremoved);
  }
  return newfreqs.head(size - nremoved);
}

writeHessianFile()

void eonc::Hessian::writeHessianFile ( bool on )

inlinenoexcept

Whether a finished Hessian goes to hessian.dat (on by default).

A caller that takes many Hessians, one per instanton bead, turns it off.

Definition at line 61 of file Hessian.h.

{ writeFile = on; }

Member Data Documentation

atoms

VectorXi eonc::Hessian::atoms

private

Definition at line 72 of file Hessian.h.

freqs

VectorXd eonc::Hessian::freqs

private

Definition at line 68 of file Hessian.h.

hessian

MatrixXd eonc::Hessian::hessian

private

Definition at line 67 of file Hessian.h.

log

eonc::log::Scoped eonc::Hessian::log

private

Definition at line 78 of file Hessian.h.

matter

Matter* eonc::Hessian::matter

private

Definition at line 64 of file Hessian.h.

modes

MatrixXd eonc::Hessian::modes

private

Definition at line 69 of file Hessian.h.

parameters

const Parameters& eonc::Hessian::parameters

private

Definition at line 65 of file Hessian.h.

writeFile

bool eonc::Hessian::writeFile = true

private

Definition at line 70 of file Hessian.h.

---

The documentation for this class was generated from the following files:

• /home/runner/work/eOn/eOn/include/eon/Hessian.h
• /home/runner/work/eOn/eOn/client/Hessian.cpp