eonc::Davidson Class Reference

Davidson method for the lowest Hessian curvature mode (min-mode). More...

#include <Davidson.h>

Inheritance diagram for eonc::Davidson:

Public Member Functions
	Davidson (std::shared_ptr< Matter > matter, const Parameters &params, std::shared_ptr< Potential > pot)
	~Davidson ()=default
void	compute (std::shared_ptr< Matter > matter, AtomMatrix initialDirection) override
void	compute (std::shared_ptr< Matter > matter, AtomMatrix initialDirection, const VectorXi &mobileAtoms)
double	getEigenvalue () override
AtomMatrix	getEigenvector () override
Public Member Functions inherited from eonc::LowestEigenmode
	LowestEigenmode (std::shared_ptr< Potential > potPassed, const Parameters &parameters)
virtual	~LowestEigenmode ()=default

Private Attributes
AtomMatrix	lowestEv
double	lowestEw
eonc::log::Scoped	log

Additional Inherited Members
Public Attributes inherited from eonc::LowestEigenmode
long	totalForceCalls {0}
double	statsTorque {0.0}
double	statsCurvature {0.0}
double	statsAngle {0.0}
long	statsRotations {0}
long	totalIterations {0}
Static Public Attributes inherited from eonc::LowestEigenmode
static const char	MINMODE_DIMER [] = "dimer"
static const char	MINMODE_GPRDIMER [] = "gprdimer"
static const char	MINMODE_LANCZOS [] = "lanczos"
static const char	MINMODE_DAVIDSON [] = "davidson"
Protected Attributes inherited from eonc::LowestEigenmode
std::shared_ptr< Potential >	pot
const Parameters &	params

Detailed Description

Davidson method for the lowest Hessian curvature mode (min-mode).

Uses the same finite-difference Hessian-vector products as Lanczos, but expands a Ritz subspace with residual correction (optionally diagonal- preconditioned) instead of dimer rotational constrained minimization. See Olsen et al., J. Chem. Phys. 121, 9776 (2004) for FD context; Davidson subspace expansion follows the classical residual-correction form.

Default Ritz space = all free atoms ([Davidson] phva_atoms = All). Optional mobile list is the PHVA active set; free/fixed stays the optimizer mask.

Definition at line 31 of file Davidson.h.

Constructor & Destructor Documentation

Davidson()

eonc::Davidson::Davidson	(	std::shared_ptr< Matter >	matter,
		const Parameters &	params,
		std::shared_ptr< Potential >	pot )

Definition at line 36 of file Davidson.cpp.

    : LowestEigenmode(pot, params) {
  lowestEv.resize(matter->numberOfAtoms(), 3);
  lowestEv.setZero();
  lowestEw = 0.0;
}

~Davidson()

eonc::Davidson::~Davidson ( )

default

Member Function Documentation

compute() [1/2]

void eonc::Davidson::compute ( std::shared_ptr< Matter > matter, AtomMatrix initialDirection )

overridevirtual

Implements eonc::LowestEigenmode.

Definition at line 44 of file Davidson.cpp.

                                                                         {
  const VectorXi mobile =
      resolveMobileAtoms(matter.get(), params.davidson_options().phva_atoms);
  compute(std::move(matter), std::move(direction), mobile);
}

compute() [2/2]

void eonc::Davidson::compute	(	std::shared_ptr< Matter >	matter,
		AtomMatrix	initialDirection,
		const VectorXi &	mobileAtoms )

Definition at line 50 of file Davidson.cpp.

                                                 {
  totalForceCalls = 0;
  statsRotations = 0;
  lowestEv.resize(matter->numberOfAtoms(), 3);
  lowestEv.setZero();
 
  const VectorXi mobile = resolveMobileAtoms(matter.get(), mobileIn);
  const int size = 3 * static_cast<int>(mobile.size());
  if (size == 0) {
    lowestEw = 0.0;
    return;
  }
 
  const long maxIter = std::max(1L, params.davidson_options().max_iterations);
  const double tol = params.davidson_options().tolerance;
  const double dr = params.main_options().finiteDifference;
  const FdScheme requested = parseFdScheme(params.hessian_options().fd_scheme);
  // Fourth-order is the only scheme that changes the matrix-free product.
  // one_sided and central keep the forward difference used by min-mode.
  const FdScheme hvpScheme =
      requested == FdScheme::Fourth ? FdScheme::Fourth : FdScheme::OneSided;
  const bool useDiagPrec = params.davidson_options().diagonal_preconditioner;
 
  MatrixXd V(size, maxIter);
  MatrixXd HV(size, maxIter);
  V.setZero();
  HV.setZero();
 
  VectorXd r = packMobileRows(direction, mobile);
  double beta = r.norm();
  if (beta < eonc::safemath::eps) {
    lowestEw = 0.0;
    return;
  }
  r /= beta;
 
  auto tmpMatter = std::make_unique<Matter>(*matter);
  const long forceCallsStart = tmpMatter->getForceCalls();
  const AtomMatrix pos0 = matter->getPositions();
  const VectorXd force0 = mobileForces(tmpMatter.get(), mobile);
 
  auto applyH = [&](const VectorXd &v) -> VectorXd {
    auto at = [&](double scale) -> VectorXd {
      AtomMatrix pos = pos0;
      unpackMobileRows(packMobileRows(pos0, mobile) + (scale * dr) * v, mobile,
                       pos);
      tmpMatter->setPositions(pos);
      return mobileForces(tmpMatter.get(), mobile);
    };
    return fdHessianVector(hvpScheme, dr, force0, at);
  };
 
  VectorXd diagH = VectorXd::Ones(size);
 
  V.col(0) = r;
  HV.col(0) = applyH(V.col(0));
  if (useDiagPrec) {
    for (int k = 0; k < size; ++k) {
      const double vk = V(k, 0);
      if (std::fabs(vk) > 1e-8) {
        diagH(k) = std::max(std::fabs(HV(k, 0) / vk), 1e-3);
      }
    }
  }
 
  double ew = 0.0, ewOld = 0.0;
  VectorXd evEst = V.col(0);
  VectorXd evOldEst = evEst;
  int subspace = 1;
 
  for (int iter = 0; iter < maxIter; ++iter) {
    statsRotations = iter;
    statsAngle = 0.0;
 
    MatrixXd G = V.leftCols(subspace).transpose() * HV.leftCols(subspace);
    G = 0.5 * (G + G.transpose());
 
    Eigen::SelfAdjointEigenSolver<MatrixXd> es(G);
    ew = es.eigenvalues()(0);
    VectorXd y = es.eigenvectors().col(0);
    evEst = V.leftCols(subspace) * y;
    eonc::safemath::safe_normalize_inplace(evEst);
 
    VectorXd Hx = HV.leftCols(subspace) * y;
    VectorXd resid = Hx - ew * evEst;
    const double residNorm = resid.norm();
    const double ewAbsRelErr =
        (iter == 0) ? 1.0
                    : eonc::safemath::safe_div(std::fabs(ew - ewOld),
                                               std::fabs(ewOld), 1.0);
    ewOld = ew;
    statsTorque = std::max(ewAbsRelErr, residNorm / (std::fabs(ew) + 1e-12));
    statsAngle = eonc::safemath::safe_acos(std::fabs(evEst.dot(evOldEst))) *
                 (180 / eonc::helpers::pi);
    evOldEst = evEst;
 
    QUILL_LOG_INFO(log,
                   "[Davidson] ew={:10.6f} rel_err={:10.6f} |r|={:10.6f} "
                   "angle={:7.3f} dim={:3d} iter={:3d} n_mobile={}",
                   ew, ewAbsRelErr, residNorm, statsAngle, subspace, iter,
                   mobile.size());
 
    if (ewAbsRelErr < tol && residNorm < tol * (std::fabs(ew) + 1.0)) {
      QUILL_LOG_INFO(log, "[Davidson] Tolerance reached: {}", tol);
      break;
    }
    if (subspace >= maxIter) {
      QUILL_LOG_ERROR(log, "[Davidson] Max subspace dimension");
      break;
    }
 
    VectorXd t(size);
    if (useDiagPrec) {
      for (int k = 0; k < size; ++k) {
        const double denom = diagH(k) - ew;
        t(k) = resid(k) / (std::fabs(denom) > 1e-12 ? denom : 1e-12);
      }
    } else {
      t = resid;
    }
 
    for (int c = 0; c < subspace; ++c) {
      t -= V.col(c).dot(t) * V.col(c);
    }
    const double tnorm = t.norm();
    if (tnorm < 1e-14) {
      QUILL_LOG_ERROR(log, "[Davidson] Linear dependence in residual");
      break;
    }
    t /= tnorm;
 
    V.col(subspace) = t;
    HV.col(subspace) = applyH(t);
    if (useDiagPrec) {
      for (int k = 0; k < size; ++k) {
        const double vk = t(k);
        if (std::fabs(vk) > 1e-8) {
          const double d = std::fabs(HV(k, subspace) / vk);
          diagH(k) = std::max(diagH(k), std::max(d, 1e-3));
        }
      }
    }
    ++subspace;
 
    if (iter >= maxIter - 1) {
      QUILL_LOG_ERROR(log, "[Davidson] Max iterations");
      break;
    }
  }
 
  lowestEw = ew;
  totalForceCalls = tmpMatter->getForceCalls() - forceCallsStart;
  lowestEv.setZero();
  if (evEst.size() == size) {
    unpackMobileRows(evEst, mobile, lowestEv);
  }
}

getEigenvalue()

double eonc::Davidson::getEigenvalue ( )

overridevirtual

Implements eonc::LowestEigenmode.

Definition at line 209 of file Davidson.cpp.

{ return lowestEw; }

getEigenvector()

AtomMatrix eonc::Davidson::getEigenvector ( )

overridevirtual

Implements eonc::LowestEigenmode.

Definition at line 211 of file Davidson.cpp.

{ return lowestEv; }

Member Data Documentation

log

eonc::log::Scoped eonc::Davidson::log

private

Definition at line 47 of file Davidson.h.

lowestEv

AtomMatrix eonc::Davidson::lowestEv

private

Definition at line 45 of file Davidson.h.

lowestEw

double eonc::Davidson::lowestEw

private

Definition at line 46 of file Davidson.h.

---

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

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