Dimer Class Reference

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

#include <Dimer.h>

Inheritance diagram for 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)
double	getEigenvalue ()
AtomMatrix	getEigenvector ()
Public Member Functions inherited from eonc::LowestEigenmode
	LowestEigenmode (std::shared_ptr< Potential > potPassed, const Parameters &parameters)
	~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()

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

Definition at line 23 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 Dimer::calcRotationalForceReturnCurvature ( AtomMatrix & rotationalForce )

private

Compute rotational force and return curvature along the dimer.

Definition at line 150 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);
    rotationRemove(matterCenter, matterDimer);
    posDimer = matterDimer->getPositions();
    direction = 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);
      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);
      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 = makeOrthogonal(forceA, direction);
  forceB = 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 Dimer::compute	(	std::shared_ptr< Matter >	matter,
		AtomMatrix	initialDirection )

Definition at line 43 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);
    *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();
    assert(std::isnormal(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;
    }
 
    double rotForce1 = matDot(rotationalForce, rotationalPlane);
    rotate(params.dimer_options.rotation_angle);
 
    if (!doneRotating) {
      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);
  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 Dimer::determineRotationalPlane ( const AtomMatrix & rotationalForce, AtomMatrix & rotationalForceOld, const AtomMatrix & rotationalPlaneOld, double & lengthRotationalForceOld )

private

Determine rotational plane via conjugate gradient.

Definition at line 232 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 = makeOrthogonal(rotationalPlane, direction);
  eonc::safemath::safe_normalize_inplace(rotationalPlane);
 
  rotationalForceOld = rotationalForce;
}

getEigenvalue()

double Dimer::getEigenvalue ( )

nodiscard

Definition at line 146 of file Dimer.cpp.

{ return eigenvalue; }

getEigenvector()

AtomMatrix Dimer::getEigenvector ( )

nodiscard

Definition at line 148 of file Dimer.cpp.

{ return direction; }

rotate()

void Dimer::rotate ( double rotationAngle )

private

Rotate the dimer by the given angle (radians).

Definition at line 259 of file Dimer.cpp.

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

Member Data Documentation

direction

AtomMatrix eonc::Dimer::direction

private

Definition at line 36 of file Dimer.h.

eigenvalue

double eonc::Dimer::eigenvalue {0.0}

private

Definition at line 38 of file Dimer.h.

{0.0};

log

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

private

Definition at line 33 of file Dimer.h.

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

matterCenter

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

private

Definition at line 34 of file Dimer.h.

matterDimer

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

private

Definition at line 35 of file Dimer.h.

nAtoms

long eonc::Dimer::nAtoms {0}

private

Definition at line 39 of file Dimer.h.

{0};

rotationalPlane

AtomMatrix eonc::Dimer::rotationalPlane

private

Definition at line 37 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