eonc::Dimer Class Reference

Classic dimer method to find the lowest curvature mode. More...

#include <Dimer.h>

Inheritance diagram for eonc::Dimer:

Public Member Functions
	Dimer (std::shared_ptr< Matter > matter, const Parameters &params, std::shared_ptr< Potential > pot)
	~Dimer ()=default
void	compute (std::shared_ptr< Matter > matter, AtomMatrix initialDirection) override
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 Member Functions
void	determineRotationalPlane (const AtomMatrix &rotationalForce, AtomMatrix &rotationalForceOld, const AtomMatrix &rotationalPlaneOld, double &lengthRotationalForceOld)
	Determine rotational plane via conjugate gradient.
void	rotate (double rotationAngle)
	Rotate the dimer by the given angle (radians).
double	calcRotationalForceReturnCurvature (AtomMatrix &rotationalForce)
	Compute rotational force and return curvature along the dimer.

Private Attributes
eonc::log::FileScoped	log {"dimer", "dimer.log"}
std::shared_ptr< Matter >	matterCenter
std::shared_ptr< Matter >	matterDimer
AtomMatrix	direction
AtomMatrix	rotationalPlane
double	eigenvalue {0.0}
long	nAtoms {0}

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

Classic dimer method to find the lowest curvature mode.

Uses finite-difference rotation to converge on the minimum eigenmode.

Definition at line 22 of file Dimer.h.

Constructor & Destructor Documentation

Dimer()

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

Definition at line 44 of file Dimer.cpp.

    : LowestEigenmode(pot, params) {
  // Give matterDimer its own potential for parallel force evaluation
  auto dimerPot =
      (pot->needsPerImageInstance() && params.main_options().parallel)
          ? eonc::helpers::makePotential(params)
          : pot;
  matterCenter = std::make_shared<Matter>(pot, params);
  matterDimer = std::make_shared<Matter>(dimerPot, params);
  *matterCenter = *matter;
  *matterDimer = *matter;
  nAtoms = matter->numberOfAtoms();
 
  direction.resize(nAtoms, 3);
  rotationalPlane.resize(nAtoms, 3);
  direction.setZero();
  rotationalPlane.setZero();
  totalForceCalls = 0;
}

~Dimer()

eonc::Dimer::~Dimer ( )

default

Member Function Documentation

calcRotationalForceReturnCurvature()

double eonc::Dimer::calcRotationalForceReturnCurvature ( AtomMatrix & rotationalForce )

private

Compute rotational force and return curvature along the dimer.

Definition at line 185 of file Dimer.cpp.

                                                                            {
  AtomMatrix posCenter = matterCenter->getPositions();
 
  // Displace to get dimer configuration A
  AtomMatrix posDimer =
      posCenter + direction * params.main_options().finiteDifference;
 
  // Optional rotation removal (Melander, Laasonen, Jonsson, JCTC 2015)
  if (params.dimer_options().remove_rotation) {
    matterDimer->setPositions(posDimer);
    eonc::geometry::rotationRemove(matterCenter, matterDimer);
    posDimer = matterDimer->getPositions();
    direction = matterCenter->pbc(posDimer - posCenter);
    eonc::safemath::safe_normalize_inplace(direction);
    posDimer = posCenter + direction * params.main_options().finiteDifference;
  }
 
  // Obtain forces for dimer and center
  matterDimer->setPositions(posDimer);
  AtomMatrix forceA, forceCenter;
  if (pot->supportsBatchEvaluation()) {
    // Only batch systems that actually need recomputation.
    bool centerDirty = matterCenter->needsForceUpdate();
    bool dimerDirty = matterDimer->needsForceUpdate();
 
    if (centerDirty && dimerDirty) {
      // Both need eval -- batch together
      auto nrs0 = matterCenter->getAtomicNrs();
      auto nrs1 = matterDimer->getAtomicNrs();
      auto box0 = matterCenter->getCell();
      auto box1 = matterDimer->getCell();
      const double *posVec[] = {matterCenter->getPositions().data(),
                                posDimer.data()};
      const int *nrsVec[] = {nrs0.data(), nrs1.data()};
      double *frcVec[] = {matterCenter->forcesData(),
                          matterDimer->forcesData()};
      double energies[2], vars[2];
      const double *boxVec[] = {box0.data(), box1.data()};
      pot->forceBatch(2, nAtoms, posVec, nrsVec, frcVec, energies, vars,
                      boxVec);
      rejectNonFiniteBatch(nAtoms, frcVec[0], energies[0]);
      rejectNonFiniteBatch(nAtoms, frcVec[1], energies[1]);
      matterCenter->setComputedPotential(energies[0], vars[0]);
      matterDimer->setComputedPotential(energies[1], vars[1]);
    } else if (dimerDirty) {
      // Only dimer moved -- eval just dimer, center is cached
      auto nrs = matterDimer->getAtomicNrs();
      auto box = matterDimer->getCell();
      const double *posVec[] = {posDimer.data()};
      const int *nrsVec[] = {nrs.data()};
      double *frcVec[] = {matterDimer->forcesData()};
      double energies[1], vars[1];
      const double *boxVec[] = {box.data()};
      pot->forceBatch(1, nAtoms, posVec, nrsVec, frcVec, energies, vars,
                      boxVec);
      rejectNonFiniteBatch(nAtoms, frcVec[0], energies[0]);
      matterDimer->setComputedPotential(energies[0], vars[0]);
    } else if (centerDirty) {
      // Only center moved (rare)
      matterCenter->getForces(); // through computePotential
    }
    // else: both cached, nothing to do
    forceCenter = matterCenter->getForces();
    forceA = matterDimer->getForces();
  } else {
    forceA = matterDimer->getForces();
    forceCenter = matterCenter->getForces();
  }
  AtomMatrix forceB = 2.0 * forceCenter - forceA;
 
  double projA = matDot(direction, forceA);
  double projB = matDot(direction, forceB);
 
  // Remove force component parallel to dimer
  forceA = helpers::makeOrthogonal(forceA, direction);
  forceB = helpers::makeOrthogonal(forceB, direction);
 
  // Rotational force = orthogonal force difference
  rotationalForce =
      (forceA - forceB) / (2.0 * params.main_options().finiteDifference);
 
  // Curvature along the dimer
  return (projB - projA) / (2.0 * params.main_options().finiteDifference);
}

compute()

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

overridevirtual

Implements eonc::LowestEigenmode.

Definition at line 65 of file Dimer.cpp.

                                                 {
  *matterCenter = *matter;
 
  eonc::safemath::safe_normalize_inplace(initialDirection);
  direction = initialDirection;
 
  // Optional: LOR / Lanczos / Davidson (enum dispatch; classical falls
  // through).
  if (auto alt = runAlternativeRotation(params.dimer_options().rotation_backend,
                                        matter, params, pot, direction,
                                        static_cast<quill::Logger *>(log))) {
    eigenvalue = alt->eigenvalue;
    direction = alt->eigenvector;
    totalForceCalls += alt->forceCalls;
    statsRotations = alt->rotations;
    eonc::safemath::safe_normalize_inplace(direction);
    for (long i = 0; i < nAtoms; ++i) {
      if (matterCenter->getFixed(i)) {
        direction.row(i).setZero();
      }
    }
    *matterCenter = *matter;
    return;
  }
 
  long rotations = 0;
  long forceCallsCenter = matterCenter->getForceCalls();
  long forceCallsDimer = matterDimer->getForceCalls();
  double curvature = 0.0;
  double rotationAngle = 0.0;
  double torque = 0.0;
 
  AtomMatrix rotationalForce(nAtoms, 3);
  AtomMatrix rotationalForceOld(nAtoms, 3);
  AtomMatrix rotationalPlaneOld(nAtoms, 3);
  rotationalForce.setZero();
  rotationalForceOld.setZero();
  rotationalPlaneOld.setZero();
 
  statsAngle = 0;
  double lengthRotationalForceOld = 0.0;
 
  // Two force calls per rotation iteration
  bool doneRotating = false;
  while (!doneRotating) {
    curvature = calcRotationalForceReturnCurvature(rotationalForce);
 
    determineRotationalPlane(rotationalForce, rotationalForceOld,
                             rotationalPlaneOld, lengthRotationalForceOld);
 
    torque = rotationalForce.norm();
    // torque == 0 is a valid aligned dimer; isnormal(0) is false.
    assert(std::isfinite(torque));
 
    // Convergence: stop if torque is below threshold or max rotations reached
    if ((torque > params.dimer_options().torque_max &&
         rotations >= params.dimer_options().rotations_max) ||
        (torque < params.dimer_options().torque_max &&
         torque >= params.dimer_options().torque_min &&
         rotations >= params.dimer_options().rotations_min) ||
        (torque < params.dimer_options().torque_min)) {
      doneRotating = true;
    }
 
    // A probe after acceptance would return rotation_angle off the
    // direction whose curvature was just measured.
    if (!doneRotating) {
      double rotForce1 = matDot(rotationalForce, rotationalPlane);
      rotate(params.dimer_options().rotation_angle);
 
      curvature = calcRotationalForceReturnCurvature(rotationalForce);
      double rotForce2 = matDot(rotationalForce, rotationalPlane);
 
      double rotForceChange =
          (rotForce1 - rotForce2) / params.dimer_options().rotation_angle;
      double forceDimer = (rotForce1 + rotForce2) / 2.0;
 
      rotationAngle = eonc::safemath::safe_atan_ratio(2.0 * forceDimer,
                                                      rotForceChange, 0.0) /
                          2.0 -
                      params.dimer_options().rotation_angle / 2.0;
 
      if (rotForceChange < 0) {
        rotationAngle += eonc::helpers::pi / 2.0;
      }
 
      rotate(rotationAngle);
      rotationalPlaneOld = rotationalPlane;
      rotations++;
    }
    QUILL_LOG_DEBUG(log,
                    "[DimerRot]   -----   ---------   ----------------   "
                    "---------  {:9.3e}  {:9.3e}  {:9.3e}   ---------\n",
                    curvature, torque,
                    rotationAngle * (180.0 / eonc::helpers::pi));
  }
 
  statsTorque = torque;
  statsCurvature = curvature;
  eonc::safemath::safe_normalize_inplace(direction);
  for (long i = 0; i < nAtoms; ++i) {
    if (matterCenter->getFixed(i)) {
      direction.row(i).setZero();
    }
  }
  statsAngle = eonc::safemath::safe_acos(matDot(direction, initialDirection));
  statsAngle *= (180.0 / eonc::helpers::pi);
  statsRotations = rotations;
  eigenvalue = curvature;
 
  forceCallsCenter = matterCenter->getForceCalls() - forceCallsCenter;
  forceCallsDimer = matterDimer->getForceCalls() - forceCallsDimer;
  totalForceCalls += forceCallsCenter + forceCallsDimer;
}

determineRotationalPlane()

void eonc::Dimer::determineRotationalPlane ( const AtomMatrix & rotationalForce, AtomMatrix & rotationalForceOld, const AtomMatrix & rotationalPlaneOld, double & lengthRotationalForceOld )

private

Determine rotational plane via conjugate gradient.

Definition at line 270 of file Dimer.cpp.

                                                                       {
  double gamma = 0.0;
  double a = std::abs(matDot(rotationalForce, rotationalForceOld));
  double b = rotationalForceOld.squaredNorm();
  if (a < 0.5 * b) {
    // Polak-Ribiere conjugate gradient direction
    gamma = (rotationalForce.array() *
             (rotationalForce - rotationalForceOld).array())
                .sum() /
            b;
  }
 
  // New rotational plane from current force and previous plane
  rotationalPlane =
      rotationalForce + rotationalPlaneOld * lengthRotationalForceOld * gamma;
 
  // Orthogonalize to dimer direction and normalize
  lengthRotationalForceOld = rotationalPlane.norm();
  rotationalPlane = helpers::makeOrthogonal(rotationalPlane, direction);
  eonc::safemath::safe_normalize_inplace(rotationalPlane);
 
  rotationalForceOld = rotationalForce;
}

getEigenvalue()

double eonc::Dimer::getEigenvalue ( )

nodiscardoverridevirtual

Implements eonc::LowestEigenmode.

Definition at line 181 of file Dimer.cpp.

{ return eigenvalue; }

getEigenvector()

AtomMatrix eonc::Dimer::getEigenvector ( )

nodiscardoverridevirtual

Implements eonc::LowestEigenmode.

Definition at line 183 of file Dimer.cpp.

{ return direction; }

rotate()

void eonc::Dimer::rotate ( double rotationAngle )

private

Rotate the dimer by the given angle (radians).

Definition at line 297 of file Dimer.cpp.

                                       {
  statsAngle += rotationAngle;
 
  double cosA = std::cos(rotationAngle);
  double sinA = std::sin(rotationAngle);
 
  AtomMatrix newDirection = direction * cosA + rotationalPlane * sinA;
  AtomMatrix newPlane = rotationalPlane * cosA - direction * sinA;
  direction = std::move(newDirection);
  rotationalPlane = std::move(newPlane);
 
  eonc::safemath::safe_normalize_inplace(direction);
  eonc::safemath::safe_normalize_inplace(rotationalPlane);
 
  // Remove component from rotationalPlane parallel to direction
  rotationalPlane = helpers::makeOrthogonal(rotationalPlane, direction);
  eonc::safemath::safe_normalize_inplace(rotationalPlane);
}

Member Data Documentation

direction

AtomMatrix eonc::Dimer::direction

private

Definition at line 37 of file Dimer.h.

eigenvalue

double eonc::Dimer::eigenvalue {0.0}

private

Definition at line 39 of file Dimer.h.

{0.0};

log

eonc::log::FileScoped eonc::Dimer::log {"dimer", "dimer.log"}

private

Definition at line 34 of file Dimer.h.

{"dimer", "dimer.log"};

matterCenter

std::shared_ptr<Matter> eonc::Dimer::matterCenter

private

Definition at line 35 of file Dimer.h.

matterDimer

std::shared_ptr<Matter> eonc::Dimer::matterDimer

private

Definition at line 36 of file Dimer.h.

nAtoms

long eonc::Dimer::nAtoms {0}

private

Definition at line 40 of file Dimer.h.

{0};

rotationalPlane

AtomMatrix eonc::Dimer::rotationalPlane

private

Definition at line 38 of file Dimer.h.

---

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

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