eonc::Lanczos Class Reference

#include <Lanczos.h>

Inheritance diagram for eonc::Lanczos:

Public Member Functions
	Lanczos (std::shared_ptr< Matter > matter, const Parameters &params, std::shared_ptr< Potential > pot)
	~Lanczos ()=default
void	compute (std::shared_ptr< Matter > matter, AtomMatrix initialDirection) override
void	compute (std::shared_ptr< Matter > matter, AtomMatrix initialDirection, const VectorXi &mobileAtoms)
	Same as compute(matter, dir) but with an explicit mobile atom list (intersected with free flags).
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 {0.0}
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

Definition at line 26 of file Lanczos.h.

Constructor & Destructor Documentation

Lanczos()

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

Definition at line 36 of file Lanczos.cpp.

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

~Lanczos()

eonc::Lanczos::~Lanczos ( )

default

Member Function Documentation

compute() [1/2]

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

overridevirtual

Implements eonc::LowestEigenmode.

Definition at line 44 of file Lanczos.cpp.

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

compute() [2/2]

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

Same as compute(matter, dir) but with an explicit mobile atom list (intersected with free flags).

Eigenvector rows outside the list are 0.

Definition at line 50 of file Lanczos.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 maxIters = std::max(1L, params.lanczos_options().max_iterations);
  MatrixXd T(size, maxIters), Q(size, maxIters);
  T.setZero();
  VectorXd u(size), r = packMobileRows(direction, mobile);
 
  double alpha, beta = r.norm();
  if (beta < eonc::safemath::eps) {
    lowestEw = 0.0;
    return;
  }
  double ew = 0, ewOld = 0, ewAbsRelErr;
  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;
  VectorXd evEst, evT, evOldEst;
 
  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);
  };
 
  for (int i = 0; i < size; i++) {
    statsRotations = i;
    Q.col(i) = r / beta;
 
    u = applyH(Q.col(i));
 
    if (i == 0) {
      r = u;
    } else {
      r = u - beta * Q.col(i - 1);
    }
    alpha = Q.col(i).dot(r);
    r = r - alpha * Q.col(i);
 
    T(i, i) = alpha;
    if (i > 0) {
      T(i - 1, i) = beta;
      T(i, i - 1) = beta;
    }
 
    beta = r.norm();
    // A vanished residual closes this Krylov block. Take the Ritz pair from
    // the current T; breaking first keeps the previous subspace.
    const bool krylovClosed = beta <= 1e-10 * std::fabs(alpha);
 
    if (i >= 1) {
      Eigen::SelfAdjointEigenSolver<MatrixXd> es(T.block(0, 0, i + 1, i + 1));
      ew = es.eigenvalues()(0);
      evT = es.eigenvectors().col(0);
      ewAbsRelErr = eonc::safemath::safe_div(std::fabs(ew - ewOld),
                                             std::fabs(ewOld), 1.0);
      ewOld = ew;
 
      evEst = Q.block(0, 0, size, i + 1) * evT;
      eonc::safemath::safe_normalize_inplace(evEst);
      statsAngle = eonc::safemath::safe_acos(std::fabs(evEst.dot(evOldEst))) *
                   (180 / eonc::helpers::pi);
      statsTorque = ewAbsRelErr;
      evOldEst = evEst;
      QUILL_LOG_INFO(log,
                     "[ILanczos] {:9s} {:9s} {:10s} {:14s} {:9.4f} "
                     "{:10.6f} {:7.3f} {:5} n_mobile={}",
                     "----", "----", "----", "----", ew, ewAbsRelErr,
                     statsAngle, i, mobile.size());
      if (krylovClosed) {
        QUILL_LOG_ERROR(log, "[ILanczos] ERROR: linear dependence");
        break;
      }
      if (ewAbsRelErr < params.lanczos_options().tolerance) {
        QUILL_LOG_INFO(log, "[ILanczos] Tolerance reached: {}",
                       params.lanczos_options().tolerance);
        break;
      }
    } else {
      ew = alpha;
      ewOld = ew;
      evEst = Q.col(0);
      evOldEst = Q.col(0);
      if (krylovClosed) {
        QUILL_LOG_ERROR(log, "[ILanczos] ERROR: linear dependence");
        break;
      }
      if (lowestEw != 0.0 && params.lanczos_options().quit_early) {
        double Cprev = lowestEw;
        double Cnew = u.dot(Q.col(i));
        ewAbsRelErr = eonc::safemath::safe_div(std::fabs(Cnew - Cprev),
                                               std::fabs(Cprev), 1.0);
        if (ewAbsRelErr <= params.lanczos_options().tolerance) {
          statsAngle = 0.0;
          statsTorque = ewAbsRelErr;
          QUILL_LOG_INFO(log, "[ILanczos] Tolerance reached: {}",
                         params.lanczos_options().tolerance);
          break;
        }
      }
    }
 
    if (i >= params.lanczos_options().max_iterations - 1) {
      QUILL_LOG_ERROR(log, "[ILanczos] Max iterations");
      break;
    }
  }
 
  lowestEw = ew;
  totalForceCalls = tmpMatter->getForceCalls() - forceCallsStart;
 
  lowestEv.setZero();
  if (evEst.size() == size) {
    unpackMobileRows(evEst, mobile, lowestEv);
  }
}

getEigenvalue()

double eonc::Lanczos::getEigenvalue ( )

overridevirtual

Implements eonc::LowestEigenmode.

Definition at line 191 of file Lanczos.cpp.

{ return lowestEw; }

getEigenvector()

AtomMatrix eonc::Lanczos::getEigenvector ( )

overridevirtual

Implements eonc::LowestEigenmode.

Definition at line 193 of file Lanczos.cpp.

{ return lowestEv; }

Member Data Documentation

log

eonc::log::Scoped eonc::Lanczos::log

private

Definition at line 44 of file Lanczos.h.

lowestEv

AtomMatrix eonc::Lanczos::lowestEv

private

Definition at line 42 of file Lanczos.h.

lowestEw

double eonc::Lanczos::lowestEw {0.0}

private

Definition at line 43 of file Lanczos.h.

{0.0};

---

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

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