eonc::neb Namespace Reference

Namespaces
namespace	zoom

Classes
struct	CartesianStep
	One steepest step in Cartesian coordinates and the cell matrix. More...
struct	DNEB_Projection
	Trygubenko & Wales, JCP 120:2082, 2004. More...
struct	ExtremaResult
struct	ImageForceData
	Data for a single image needed by projection strategies. More...
struct	ImprovedTangent
	Henkelman & Jonsson, JCP 113:9978, 2000. More...
struct	JointBlock
struct	NEB_Projection
	Jonsson, Mills, Jacobsen 1998 (World Scientific). More...
class	OCINEBController
	Goswami (in prep). More...
struct	OnsagerMachlupSpring
	Onsager-Machlup action-based springs (Mandelli & Parrinello 2021). More...
struct	PlainEB
	Mills, Jonsson, Schenter, Surf. More...
struct	SimpleTangent
	Mills, Jonsson, Schenter, Surf. More...
struct	SpringResult
	Result of spring force computation for a single image. More...
struct	UniformSpring
	Uniform spring constant for all images. More...
struct	WeightedSpring
	Energy-weighted spring constants (variable per segment). More...

Typedefs
using	ProjectionStrategy
using	SpringStrategy
using	TangentStrategy = std::variant<SimpleTangent, ImprovedTangent>

Functions
AtomMatrix	computeTangent (const AtomMatrix &posDiffNext, const AtomMatrix &posDiffPrev, double energy, double energyPrev, double energyNext, bool use_old_tangent)
	Compute the tangent vector at image i using the improved tangent scheme.
AtomMatrix	forcePerp (const AtomMatrix &force, const AtomMatrix &tangent)
	Compute the perpendicular component of force relative to the tangent.
AtomMatrix	climbingImageForce (const AtomMatrix &force, const AtomMatrix &tangent, const AtomMatrix &forceDNEB)
	Compute the climbing image projected force.
AtomMatrix	computeDNEB (const AtomMatrix &forceSpring, const AtomMatrix &tangent, const AtomMatrix &forcePerp, bool useSwitching=true)
	Compute the doubly-nudged elastic band perpendicular spring force.
void	zeroTranslation (AtomMatrix &projectedForce, int nFreeAtoms, int nAtoms)
	Zero net translational force for fully free systems.
ProjectionStrategy	buildProjectionStrategy (const Parameters &params)
	Build the projection strategy from parameters.
AtomMatrix	computeDNEBComponent (const AtomMatrix &forceSpring, const AtomMatrix &tangent, const AtomMatrix &forcePerp, bool useSwitching)
	Compute the DNEB force component for a given image.
ExtremaResult	findSplineExtrema (const std::vector< std::shared_ptr< Matter > > &path, const std::vector< std::shared_ptr< AtomMatrix > > &tangent, long numImages)
	Find extrema along the MEP using cubic spline interpolation.
AtomMatrix	interpolatedPeakMode (const std::vector< std::shared_ptr< Matter > > &path, const std::vector< std::shared_ptr< AtomMatrix > > &tangent, long numImages, double posFraction)
	Unit tangent at a fractional image index.
void	printImageData (const std::vector< std::shared_ptr< Matter > > &path, const std::vector< std::shared_ptr< AtomMatrix > > &tangent, const std::vector< std::shared_ptr< EigenmodeStrategy > > &eigenmode_solvers, long numImages, bool estimateEigenvalues, bool writeToFile, size_t idx, eonc::log::Scoped log, double referenceEnergy=std::numeric_limits< double >::quiet_NaN())
	Print NEB image data to log and optionally to file.
std::vector< readcon::ConFrame >	pathToConFrames (const std::vector< std::shared_ptr< Matter > > &path, const std::vector< std::shared_ptr< AtomMatrix > > &tangent, const std::vector< std::shared_ptr< EigenmodeStrategy > > &eigenmode_solvers, long numImages, bool estimateEigenvalues, std::optional< size_t > bandIndex=std::nullopt, double referenceEnergy=std::numeric_limits< double >::quiet_NaN())
	Build stamped ConFrames for a NEB band (same metadata as writePathCon).
eonc::io::IoStatus	writePathCon (const std::vector< std::shared_ptr< Matter > > &path, const std::vector< std::shared_ptr< AtomMatrix > > &tangent, const std::vector< std::shared_ptr< EigenmodeStrategy > > &eigenmode_solvers, long numImages, bool estimateEigenvalues, std::string filename, std::optional< size_t > bandIndex=std::nullopt, double referenceEnergy=std::numeric_limits< double >::quiet_NaN())
	Write a NEB band as a multi-frame .con via readcon ConFrameBuilder::clone().
SpringStrategy	buildSpringStrategy (const Parameters &params, const std::vector< std::shared_ptr< Matter > > &path, long numImages, int atoms, double maxEnergy, double E_ref)
	Build the appropriate spring strategy from parameters and current path state.
TangentStrategy	buildTangentStrategy (const Parameters &params)
	Build the tangent strategy from parameters.
double	solidStateJacobian (double meanVolume, int nAtoms, double weight)
	J = (V/N)^{1/3} * N^{1/2} * weight, with V the mean endpoint volume.
bool	orientCellLowerTriangular (Matrix3d &cell, AtomMatrix &positions)
	Rotate the lab frame so the cell is lower triangular: the first lattice vector lies on x and the second lies in the xy plane.
void	orientSolidStateMatter (Matter &image)
void	interpolateSolidStateLinear (std::vector< Matter > &images)
	Replace interior images by a fractional linear interpolation of the oriented endpoints.
JointBlock	jointDisplacement (const Matter &from, const Matter &to, double jacobian)
	Displacement of to relative to from in the joint metric.
double	jointNorm (const JointBlock &block)
Matrix3d	cellNebForce (const Matrix3d &cauchy, double volume, double jacobian, const Matrix3d &externalStress)
	NEB force on the Jacobian-scaled strain.
Matrix3d	finiteDifferenceCauchyStress (const Matter &image, double strainStep)
	Central difference of the potential energy on the six lower strain components.
double	solidStateEnthalpy (const Matter &image, const Matter &reference, double pressure)
	Potential energy plus P : (h0^{-1} (h-h0)) * V0.
CartesianStep	solidStateCartesianStep (const Matter &image, const AtomMatrix &atomicForce, const Matrix3d &cellForce, double jacobian)

Typedef Documentation

ProjectionStrategy

using eonc::neb::ProjectionStrategy

Initial value:

 
std::variant<PlainEB, NEB_Projection, DNEB_Projection>

Definition at line 50 of file NEBProjection.h.

SpringStrategy

using eonc::neb::SpringStrategy

Initial value:

 
std::variant<UniformSpring, WeightedSpring, OnsagerMachlupSpring>

Definition at line 58 of file NEBSpringForce.h.

TangentStrategy

using eonc::neb::TangentStrategy = std::variant<SimpleTangent, ImprovedTangent>

Definition at line 38 of file NEBTangent.h.

Function Documentation

buildProjectionStrategy()

ProjectionStrategy eonc::neb::buildProjectionStrategy

(

const Parameters &

params

)

Build the projection strategy from parameters.

DNEB is incompatible with OM and weighted springs; falls back to NEB_Projection.

Definition at line 33 of file NEBProjection.cpp.

                                                                     {
  bool omActive = params.neb_options().spring.om.enabled;
  bool weightedActive = params.neb_options().spring.weighting.enabled;
 
  if (params.neb_options().spring.use_elastic_band && !omActive &&
      !weightedActive) {
    return PlainEB{};
  }
  if (params.neb_options().spring.doubly_nudged && !omActive &&
      !weightedActive) {
    return DNEB_Projection{params.neb_options().spring.use_switching};
  }
  return NEB_Projection{};
}

buildSpringStrategy()

SpringStrategy eonc::neb::buildSpringStrategy	(	const Parameters &	params,
		const std::vector< std::shared_ptr< Matter > > &	path,
		long	numImages,
		int	atoms,
		double	maxEnergy,
		double	E_ref )

Build the appropriate spring strategy from parameters and current path state.

Definition at line 71 of file NEBSpringForce.cpp.

                                                                               {
 
  if (params.neb_options().spring.om.enabled) {
    double base_k = params.neb_options().spring.constant;
 
    if (params.neb_options().spring.om.optimize_k) {
      double avgPotForce = 0.0;
      double avgPathCurvature = 0.0;
      int count = 0;
      for (long j = 1; j <= numImages; j++) {
        avgPotForce += path[j]->getForces().norm();
        AtomMatrix next = path[j + 1]->getPositions();
        AtomMatrix prev = path[j - 1]->getPositions();
        AtomMatrix curr = path[j]->getPositions();
        AtomMatrix curvVec = path[j]->pbc(next + prev - 2.0 * curr);
        avgPathCurvature += curvVec.norm();
        count++;
      }
      if (count > 0 && avgPathCurvature > 1e-6) {
        double scale = params.neb_options().spring.om.k_scale;
        base_k = scale * (avgPotForce / avgPathCurvature);
        base_k = std::max(base_k, params.neb_options().spring.om.k_min);
        base_k = std::min(base_k, params.neb_options().spring.om.k_max);
      }
    }
 
    // Pre-calculate L vectors for all images
    std::vector<AtomMatrix> L_vecs(numImages + 2);
    for (long j = 0; j <= numImages + 1; j++) {
      L_vecs[j].resize(atoms, 3);
      if (j == 0 || j == numImages + 1) {
        L_vecs[j].setZero();
      } else {
        const AtomMatrix &forces = path[j]->getForces();
        double alpha_k = eonc::safemath::safe_recip(2.0 * base_k, 0.0);
        for (int k = 0; k < atoms; k++) {
          L_vecs[j].row(k) = alpha_k * forces.row(k);
        }
      }
    }
 
    return OnsagerMachlupSpring{base_k, std::move(L_vecs)};
 
  } else if (params.neb_options().spring.weighting.enabled) {
    double k_l = params.neb_options().spring.weighting.k_min;
    double k_u = params.neb_options().spring.weighting.k_max;
    std::vector<double> springConstants(numImages + 2, k_l);
 
    double energyRange = maxEnergy - E_ref;
    if (energyRange < 1e-10) {
      std::fill(springConstants.begin(), springConstants.end(), k_l);
    } else {
      for (int idx = 1; idx <= numImages + 1; idx++) {
        double Ei = std::max(path[idx]->getPotentialEnergy(),
                             path[idx - 1]->getPotentialEnergy());
        if (Ei > E_ref) {
          double alpha_i = (maxEnergy - Ei) / energyRange;
          alpha_i = std::max(0.0, std::min(1.0, alpha_i));
          springConstants[idx - 1] = (1.0 - alpha_i) * k_u + alpha_i * k_l;
        } else {
          springConstants[idx - 1] = k_l;
        }
      }
    }
 
    return WeightedSpring{std::move(springConstants)};
 
  } else {
    return UniformSpring{params.neb_options().spring.constant};
  }
}

buildTangentStrategy()

TangentStrategy eonc::neb::buildTangentStrategy

(

const Parameters &

params

)

Build the tangent strategy from parameters.

Definition at line 72 of file NEBTangent.cpp.

                                                               {
  if (params.neb_options().climbing_image.use_old_tangent) {
    return SimpleTangent{};
  }
  return ImprovedTangent{};
}

cellNebForce()

Matrix3d eonc::neb::cellNebForce	(	const Matrix3d &	cauchy,
		double	volume,
		double	jacobian,
		const Matrix3d &	externalStress )

NEB force on the Jacobian-scaled strain.

Cauchy stress uses sigma = (1/V) dE/dε for h <- h (I+ε) at fixed fractional coordinates. External stress is added before the -V factor (positive hydrostatic pressure pushes the cell inward).

Definition at line 146 of file SolidStateNEB.cpp.

                                                      {
  if (!(volume > 0.0) || !(jacobian > 0.0)) {
    throw std::invalid_argument(
        "solid-state cell force needs a positive volume and Jacobian");
  }
  Matrix3d force = -(volume / jacobian) * (cauchy + externalStress);
  zeroUpper(force);
  return force;
}

climbingImageForce()

AtomMatrix eonc::neb::climbingImageForce	(	const AtomMatrix &	force,
		const AtomMatrix &	tangent,
		const AtomMatrix &	forceDNEB )

Compute the climbing image projected force.

F_CI = F - 2*(F.t)*t + forceDNEB

Definition at line 73 of file NEBForceProjection.cpp.

                                                           {
  return force - 2.0 * matDot(force, tangent) * tangent + forceDNEB;
}

computeDNEB()

AtomMatrix eonc::neb::computeDNEB	(	const AtomMatrix &	forceSpring,
		const AtomMatrix &	tangent,
		const AtomMatrix &	forcePerp,
		bool	useSwitching = true )

Compute the doubly-nudged elastic band perpendicular spring force.

useSwitching multiplies the remainder by (2/pi)*atan(|F_perp|^2 / |F_spring_perp|^2).

Definition at line 79 of file NEBForceProjection.cpp.

                                                                   {
  AtomMatrix forceSpringPerp =
      forceSpring - matDot(forceSpring, tangent) * tangent;
 
  const double forceSpringPerpNorm = forceSpringPerp.norm();
  const double forcePerpNorm = fPerp.norm();
 
  if (forceSpringPerpNorm > 1e-10 && forcePerpNorm > 1e-10) {
    AtomMatrix forcePerpNormalized = fPerp / forcePerpNorm;
    AtomMatrix dneb =
        forceSpringPerp -
        matDot(forceSpringPerp, forcePerpNormalized) * forcePerpNormalized;
 
    if (useSwitching) {
      double switching = 2.0 / eonc::helpers::pi *
                         std::atan(forcePerpNorm * forcePerpNorm /
                                   (forceSpringPerpNorm * forceSpringPerpNorm));
      dneb *= switching;
    }
    return dneb;
  }
 
  return AtomMatrix::Zero(tangent.rows(), tangent.cols());
}

computeDNEBComponent()

AtomMatrix eonc::neb::computeDNEBComponent	(	const AtomMatrix &	forceSpring,
		const AtomMatrix &	tangent,
		const AtomMatrix &	forcePerp,
		bool	useSwitching )

Compute the DNEB force component for a given image.

Returned separately so it can be added to the CI force when DNEB is active. useSwitching applies the (2/pi)*atan scale. Otherwise the perpendicular spring remainder is unscaled.

Definition at line 48 of file NEBProjection.cpp.

                                                                            {
  return eonc::neb::computeDNEB(forceSpring, tangent, fPerp, useSwitching);
}

computeTangent()

AtomMatrix eonc::neb::computeTangent	(	const AtomMatrix &	posDiffNext,
		const AtomMatrix &	posDiffPrev,
		double	energy,
		double	energyPrev,
		double	energyNext,
		bool	use_old_tangent )

Compute the tangent vector at image i using the improved tangent scheme.

Returns a normalized tangent vector.

Definition at line 22 of file NEBForceProjection.cpp.

                                                {
  AtomMatrix tang;
 
  if (use_old_tangent) {
    tang = posDiffNext;
  } else {
    // Improved tangent scheme
    if (energyNext > energy && energy > energyPrev) {
      tang = posDiffNext;
    } else if (energy > energyNext && energyPrev > energy) {
      tang = posDiffPrev;
    } else {
      // Extremum: energy-weighted combination
      double energyDiffPrev = energyPrev - energy;
      double energyDiffNext = energyNext - energy;
      double minDiffEnergy =
          std::min(std::abs(energyDiffPrev), std::abs(energyDiffNext));
      double maxDiffEnergy =
          std::max(std::abs(energyDiffPrev), std::abs(energyDiffNext));
 
      if (energyDiffPrev > energyDiffNext) {
        tang = posDiffNext * minDiffEnergy + posDiffPrev * maxDiffEnergy;
      } else {
        tang = posDiffNext * maxDiffEnergy + posDiffPrev * minDiffEnergy;
      }
    }
  }
 
  // Normalize with safety check
  double norm = tang.norm();
  if (norm > 1e-10) {
    tang /= norm;
  } else {
    // Fallback: use direction to next image
    tang = posDiffNext;
    norm = tang.norm();
    if (norm > 1e-10) {
      tang /= norm;
    }
  }
 
  return tang;
}

findSplineExtrema()

ExtremaResult eonc::neb::findSplineExtrema	(	const std::vector< std::shared_ptr< Matter > > &	path,
		const std::vector< std::shared_ptr< AtomMatrix > > &	tangent,
		long	numImages )

Find extrema along the MEP using cubic spline interpolation.

Definition at line 111 of file NEBSplineExtrema.cpp.

                                  {
 
  auto *log = eonc::log::get();
 
  // Calculate cubic parameters for each interval
  AtomMatrix tangentEndpoint;
  std::vector<double> a(numImages + 1), b(numImages + 1), c(numImages + 1),
      d(numImages + 1);
  double F1, F2, U1, U2, dist;
 
  for (long i = 0; i <= numImages; i++) {
    dist = path[i]->distanceTo(*path[i + 1]);
    if (i == 0) {
      tangentEndpoint =
          path[i]->pbc(path[1]->getPositions() - path[0]->getPositions());
      normalizeOrZero(tangentEndpoint);
      F1 = matDot(path[i]->getForces(), tangentEndpoint) * dist;
    } else {
      F1 = matDot(path[i]->getForces(), *tangent[i]) * dist;
    }
    if (i == numImages) {
      tangentEndpoint = path[i + 1]->pbc(path[numImages + 1]->getPositions() -
                                         path[numImages]->getPositions());
      normalizeOrZero(tangentEndpoint);
      F2 = matDot(path[i + 1]->getForces(), tangentEndpoint) * dist;
    } else {
      F2 = matDot(path[i + 1]->getForces(), *tangent[i + 1]) * dist;
    }
    U1 = path[i]->getPotentialEnergy();
    U2 = path[i + 1]->getPotentialEnergy();
    a[i] = U1;
    b[i] = -F1;
    c[i] = 3. * (U2 - U1) + 2. * F1 + F2;
    d[i] = -2. * (U2 - U1) - (F1 + F2);
  }
 
  ExtremaResult result;
  result.positions.resize(2 * (numImages + 1));
  result.energies.resize(2 * (numImages + 1));
  result.curvatures.resize(2 * (numImages + 1));
 
  double discriminant, f;
 
  for (long i = 0; i <= numImages; i++) {
    discriminant = c[i] * c[i] - 3.0 * b[i] * d[i];
    if (discriminant >= 0) {
      f = -1;
 
      // Quadratic case
      if ((d[i] == 0) && (c[i] != 0)) {
        f = (-b[i] / (2. * c[i]));
      }
      // Cubic case 1
      else if (d[i] != 0) {
        f = -(c[i] + std::sqrt(discriminant)) / (3. * d[i]);
      }
      if ((f >= 0) && (f <= 1)) {
        result.positions[result.numExtrema] = i + f;
        result.energies[result.numExtrema] =
            ((d[i] * f + c[i]) * f + b[i]) * f + a[i]; // Horner's method
        result.curvatures[result.numExtrema] = 6.0 * d[i] * f + 2 * c[i];
        result.numExtrema++;
      }
      // Cubic case 2. A zero cubic coefficient leaves f at the quadratic
      // root, and a zero discriminant is that same cubic root. Storing
      // either again doubles one stationary point.
      if (d[i] != 0 && discriminant != 0.0) {
        f = -(c[i] - std::sqrt(discriminant)) / (3. * d[i]);
      } else {
        f = -1;
      }
      if ((f >= 0) && (f <= 1)) {
        result.positions[result.numExtrema] = i + f;
        result.energies[result.numExtrema] =
            ((d[i] * f + c[i]) * f + b[i]) * f + a[i]; // Horner's method
        result.curvatures[result.numExtrema] = 6 * d[i] * f + 2 * c[i];
        result.numExtrema++;
      }
    }
  }
 
  QUILL_LOG_DEBUG(log, "Found {} extrema", result.numExtrema);
  QUILL_LOG_DEBUG(log, "Energy reference: {}", path[0]->getPotentialEnergy());
  for (long i = 0; i < result.numExtrema; i++) {
    QUILL_LOG_DEBUG(
        log, "extrema #{} at image position {} with energy {} and curvature {}",
        i + 1, result.positions[i],
        result.energies[i] - path[0]->getPotentialEnergy(),
        result.curvatures[i]);
  }
 
  return result;
}

finiteDifferenceCauchyStress()

Matrix3d eonc::neb::finiteDifferenceCauchyStress	(	const Matter &	image,
		double	strainStep )

Central difference of the potential energy on the six lower strain components.

Upper-triangle components stay zero.

Definition at line 157 of file SolidStateNEB.cpp.

                                                                              {
  if (!(strainStep > 0.0)) {
    throw std::invalid_argument(
        "stress finite difference needs a positive step");
  }
  const Matrix3d cell = image.getCell();
  const AtomMatrix positions = image.getPositions();
  const double volume = std::abs(cell.determinant());
  if (!(volume > 0.0)) {
    throw std::invalid_argument(
        "stress finite difference needs a nonzero cell");
  }
  Matrix3d sigma = Matrix3d::Zero();
  const int rows[6] = {0, 1, 1, 2, 2, 2};
  const int cols[6] = {0, 0, 1, 0, 1, 2};
  for (int comp = 0; comp < 6; ++comp) {
    Matrix3d strain = Matrix3d::Zero();
    strain(rows[comp], cols[comp]) = strainStep;
    const Matrix3d plus = Matrix3d::Identity() + strain;
    const Matrix3d minus = Matrix3d::Identity() - strain;
    Matter raised(image);
    raised.setCell(cell * plus);
    raised.setPositions(positions * plus);
    Matter lowered(image);
    lowered.setCell(cell * minus);
    lowered.setPositions(positions * minus);
    const double derivative =
        (raised.getPotentialEnergy() - lowered.getPotentialEnergy()) /
        (2.0 * strainStep);
    sigma(rows[comp], cols[comp]) = derivative / volume;
  }
  return sigma;
}

forcePerp()

AtomMatrix eonc::neb::forcePerp	(	const AtomMatrix &	force,
		const AtomMatrix &	tangent )

Compute the perpendicular component of force relative to the tangent.

Definition at line 69 of file NEBForceProjection.cpp.

                                                                         {
  return force - matDot(force, tangent) * tangent;
}

interpolatedPeakMode()

AtomMatrix eonc::neb::interpolatedPeakMode ( const std::vector< std::shared_ptr< Matter > > & path, const std::vector< std::shared_ptr< AtomMatrix > > & tangent, long numImages, double posFraction )

nodiscard

Unit tangent at a fractional image index.

End intervals use the geometric endpoint displacement because stored endpoint tangents stay zero.

Definition at line 208 of file NEBSplineExtrema.cpp.

                                                         {
  const int nat = path.empty() ? 0 : path[0]->numberOfAtoms();
  const auto leftIdx = static_cast<long>(std::floor(posFraction));
  if (path.size() < 2 || leftIdx < 0 || leftIdx >= numImages + 1 ||
      leftIdx + 1 >= static_cast<long>(path.size())) {
    return AtomMatrix::Zero(nat, 3);
  }
  const double f = posFraction - static_cast<double>(leftIdx);
  AtomMatrix mode =
      (1.0 - f) * storedOrEndpointTangent(path, tangent, numImages, leftIdx) +
      f * storedOrEndpointTangent(path, tangent, numImages, leftIdx + 1);
  normalizeOrZero(mode);
  return mode;
}

interpolateSolidStateLinear()

void eonc::neb::interpolateSolidStateLinear

(

std::vector< Matter > &

images

)

Replace interior images by a fractional linear interpolation of the oriented endpoints.

images includes the endpoints.

Definition at line 93 of file SolidStateNEB.cpp.

                                                            {
  const auto n = static_cast<long>(images.size());
  if (n < 3) {
    throw std::invalid_argument(
        "solid_state interpolation needs two endpoints and one image");
  }
  const Matrix3d h0 = images.front().getCell();
  const Matrix3d h1 = images.back().getCell();
  const AtomMatrix s0 = fractionalCoordinates(images.front());
  AtomMatrix ds = fractionalCoordinates(images.back()) - s0;
  for (int atom = 0; atom < ds.rows(); ++atom) {
    for (int axis = 0; axis < 3; ++axis) {
      ds(atom, axis) = wrapFrac(ds(atom, axis));
    }
  }
  const double denom = static_cast<double>(n - 1);
  for (long i = 1; i < n - 1; ++i) {
    const double t = static_cast<double>(i) / denom;
    Matrix3d h = (1.0 - t) * h0 + t * h1;
    zeroUpper(h);
    const AtomMatrix fractional = s0 + t * ds;
    images[static_cast<size_t>(i)].setCell(h);
    images[static_cast<size_t>(i)].setPositions(fractional * h);
  }
}

jointDisplacement()

JointBlock eonc::neb::jointDisplacement	(	const Matter &	from,
		const Matter &	to,
		double	jacobian )

Displacement of to relative to from in the joint metric.

Atomic rows are fractional minimum-image displacements mapped through the average cell. The cell block is the averaged right strain of the Jacobian-scaled cell difference.

Definition at line 119 of file SolidStateNEB.cpp.

                                              {
  if (!(jacobian > 0.0)) {
    throw std::invalid_argument(
        "solid-state displacement needs a positive Jacobian");
  }
  const Matrix3d hFrom = from.getCell();
  const Matrix3d hTo = to.getCell();
  AtomMatrix frac = fractionalCoordinates(to) - fractionalCoordinates(from);
  for (int atom = 0; atom < frac.rows(); ++atom) {
    for (int axis = 0; axis < 3; ++axis) {
      frac(atom, axis) = wrapFrac(frac(atom, axis));
    }
  }
  const Matrix3d average = 0.5 * (hFrom + hTo);
  JointBlock out;
  out.atomic = frac * average;
  const Matrix3d dh = jacobian * (hTo - hFrom);
  out.cell = 0.5 * (hFrom.inverse() * dh + hTo.inverse() * dh);
  zeroUpper(out.cell);
  return out;
}

jointNorm()

double eonc::neb::jointNorm

(

const JointBlock &

block

)

Definition at line 142 of file SolidStateNEB.cpp.

                                          {
  return std::hypot(block.atomic.norm(), block.cell.norm());
}

orientCellLowerTriangular()

bool eonc::neb::orientCellLowerTriangular	(	Matrix3d &	cell,
		AtomMatrix &	positions )

Rotate the lab frame so the cell is lower triangular: the first lattice vector lies on x and the second lies in the xy plane.

Fractional coordinates and the sign of the volume are unchanged. Returns false when the cell is singular or two lattice vectors are parallel.

Definition at line 49 of file SolidStateNEB.cpp.

                                                                      {
  const Vector3d a = cell.row(0);
  const Vector3d b = cell.row(1);
  const double na = a.norm();
  if (!(na > 1e-12)) {
    return false;
  }
  const Vector3d e1 = a / na;
  const Vector3d bPerp = b - b.dot(e1) * e1;
  const double nb = bPerp.norm();
  if (!(nb > 1e-12)) {
    return false;
  }
  const Vector3d e2 = bPerp / nb;
  const Vector3d e3 = e1.cross(e2);
  Matrix3d rotation;
  rotation.col(0) = e1;
  rotation.col(1) = e2;
  rotation.col(2) = e3;
  Matrix3d oriented = cell * rotation;
  zeroUpper(oriented);
  if (positions.rows() > 0) {
    positions = positions * rotation;
  }
  if (!(std::abs(oriented.determinant()) > 1e-18)) {
    return false;
  }
  cell = oriented;
  return true;
}

orientSolidStateMatter()

void eonc::neb::orientSolidStateMatter

(

Matter &

image

)

Definition at line 80 of file SolidStateNEB.cpp.

                                           {
  Matrix3d cell = image.getCell();
  AtomMatrix positions = image.getPositions();
  if (!orientCellLowerTriangular(cell, positions)) {
    throw std::invalid_argument(
        "solid_state NEB could not put the cell in lower-triangular form");
  }
  image.setCell(cell);
  if (positions.rows() > 0) {
    image.setPositions(positions);
  }
}

pathToConFrames()

std::vector< readcon::ConFrame > eonc::neb::pathToConFrames ( const std::vector< std::shared_ptr< Matter > > & path, const std::vector< std::shared_ptr< AtomMatrix > > & tangent, const std::vector< std::shared_ptr< EigenmodeStrategy > > & eigenmode_solvers, long numImages, bool estimateEigenvalues, std::optional< size_t > bandIndex = std::nullopt, double referenceEnergy = std::numeric_limits< double >::quiet_NaN() )

nodiscard

Build stamped ConFrames for a NEB band (same metadata as writePathCon).

Empty on invalid path size. Does not write to disk. referenceEnergy is the energy relative energies are measured from; NaN means path[0]. A zoomed band passes its original reactant energy, since zoom moves path[0] to the start of its window.

Definition at line 310 of file NEBSplineExtrema.cpp.

                            {
  const size_t nframes = static_cast<size_t>(numImages) + 2;
  if (path.size() < nframes) {
    return {};
  }
 
  double distTotal = 0.0;
  std::vector<eonc::io::ConFrameMetadata> metas;
  metas.reserve(nframes);
 
  // The mass-weighted arc length rides beside the Cartesian one: it is the
  // coordinate a tunnelling integral along the band runs over. A structure
  // without masses leaves it out rather than failing the write.
  double distMw = 0.0;
  bool haveMw = true;
  for (long i = 0; i <= numImages + 1; i++) {
    if (i > 0) {
      distTotal += path[i]->distanceTo(*path[i - 1]);
      if (haveMw) {
        try {
          distMw +=
              eonc::tunneling::massWeightedDistance(*path[i - 1], *path[i]);
        } catch (const std::invalid_argument &) {
          haveMw = false;
        }
      }
    }
    metas.push_back(neb_frame_metadata(path, tangent, eigenmode_solvers,
                                       numImages, estimateEigenvalues, i,
                                       distTotal, bandIndex, referenceEnergy));
    if (haveMw) {
      metas.back().scalars.push_back({"reaction_coordinate_mw", distMw});
    }
  }
  // The band's one-dimensional tunnelling estimate goes on its first frame:
  // the well frequencies, the WKB action and the splitting a two-level-system
  // screen reads. Left out when an image has no mass or an end well is flat.
  if (haveMw && !metas.empty()) {
    try {
      std::vector<std::shared_ptr<Matter>> ends(
          path.begin(), path.begin() + static_cast<std::ptrdiff_t>(nframes));
      const auto split = eonc::tunneling::bandSplitting(
          ends, referenceOrFirst(path, referenceEnergy));
      auto &head = metas.front().scalars;
      head.push_back({"hbar_omega_reactant", split.hwReactant});
      head.push_back({"hbar_omega_product", split.hwProduct});
      head.push_back({"tunnel_action", split.action});
      head.push_back({"tunnel_splitting", split.delta0});
      head.push_back({"tls_energy", split.tlsEnergy()});
      head.push_back({"tunnel_deep_wells", split.deepWells ? 1.0 : 0.0});
    } catch (const std::invalid_argument &) {
    }
  }
 
  std::vector<std::shared_ptr<Matter>> band(
      path.begin(), path.begin() + static_cast<std::ptrdiff_t>(nframes));
  return eonc::io::buildNebPathFrames(band, metas);
}

printImageData()

void eonc::neb::printImageData	(	const std::vector< std::shared_ptr< Matter > > &	path,
		const std::vector< std::shared_ptr< AtomMatrix > > &	tangent,
		const std::vector< std::shared_ptr< EigenmodeStrategy > > &	eigenmode_solvers,
		long	numImages,
		bool	estimateEigenvalues,
		bool	writeToFile,
		size_t	idx,
		eonc::log::Scoped	log,
		double	referenceEnergy )

Print NEB image data to log and optionally to file.

Definition at line 225 of file NEBSplineExtrema.cpp.

                                               {
 
  double dist, distTotal = 0;
  AtomMatrix tangentStart =
      path[0]->pbc(path[1]->getPositions() - path[0]->getPositions());
  AtomMatrix tangentEnd = path[numImages]->pbc(
      path[numImages + 1]->getPositions() - path[numImages]->getPositions());
  AtomMatrix tang;
  std::string header;
  if (estimateEigenvalues) {
    header = std::format("{:>3s} {:>12s} {:>12s} {:>12s} {:>12s}", "img",
                         "rxn_coord", "energy", "f_para", "eigval");
  } else {
    header = std::format("{:>3s} {:>12s} {:>12s} {:>12s}", "img", "rxn_coord",
                         "energy", "f_para");
  }
 
  normalizeOrZero(tangentStart);
  normalizeOrZero(tangentEnd);
 
  std::ofstream fileLogger;
  if (writeToFile) {
    std::string neb_dat_fs;
    if (idx == std::numeric_limits<size_t>::max()) {
      neb_dat_fs = "neb.dat";
    } else {
      neb_dat_fs = std::format("neb_{:03}.dat", idx);
    }
    if (fs::exists(neb_dat_fs)) {
      fs::remove(neb_dat_fs);
    }
    fileLogger.open(neb_dat_fs);
    if (fileLogger.is_open()) {
      fileLogger << header << "\n";
    }
  }
  const double energy_reactant = referenceOrFirst(path, referenceEnergy);
 
  for (long i = 0; i <= numImages + 1; i++) {
    if (i == 0) {
      tang = tangentStart;
    } else if (i == numImages + 1) {
      tang = tangentEnd;
    } else {
      tang = *tangent[i];
    }
 
    if (i > 0) {
      dist = path[i]->distanceTo(*path[i - 1]);
      distTotal += dist;
    }
 
    double relative_energy = path[i]->getPotentialEnergy() - energy_reactant;
    double parallel_force = matDot(path[i]->getForces(), tang);
 
    if (estimateEigenvalues) {
      eonc::eigenmodeCompute(*eigenmode_solvers[i], path[i], tang);
      double lowest_eigenvalue =
          eonc::eigenmodeGetEigenvalue(*eigenmode_solvers[i]);
      if (fileLogger.is_open()) {
        fileLogger << std::format(
            "{:>3} {:>12.6f} {:>12.6f} {:>12.6f} {:>12.6f}\n", i, distTotal,
            relative_energy, parallel_force, lowest_eigenvalue);
      } else {
        QUILL_LOG_DEBUG(log, "{:>3} {:>12.6f} {:>12.6f} {:>12.6f} {:>12.6f}", i,
                        distTotal, relative_energy, parallel_force,
                        lowest_eigenvalue);
      }
    } else {
      if (fileLogger.is_open()) {
        fileLogger << std::format("{:>3} {:>12.6f} {:>12.6f} {:>12.6f}\n", i,
                                  distTotal, relative_energy, parallel_force);
      } else {
        QUILL_LOG_DEBUG(log, "{:>3} {:>12.6f} {:>12.6f} {:>12.6f}", i,
                        distTotal, relative_energy, parallel_force);
      }
    }
  }
}

solidStateCartesianStep()

CartesianStep eonc::neb::solidStateCartesianStep	(	const Matter &	image,
		const AtomMatrix &	atomicForce,
		const Matrix3d &	cellForce,
		double	jacobian )

Definition at line 203 of file SolidStateNEB.cpp.

                                                       {
  if (!(jacobian > 0.0)) {
    throw std::invalid_argument("solid-state step needs a positive Jacobian");
  }
  if (atomicForce.rows() != image.numberOfAtoms()) {
    throw std::invalid_argument("solid-state step force row count mismatch");
  }
  const Matrix3d cell = image.getCell();
  Matrix3d strain = cellForce / jacobian;
  zeroUpper(strain);
  CartesianStep step;
  step.cell = cell * strain;
  zeroUpper(step.cell);
  step.positions = atomicForce + image.getPositions() * strain;
  for (long atom = 0; atom < image.numberOfAtoms(); ++atom) {
    if (image.getFixed(atom)) {
      step.positions.row(atom).setZero();
    }
  }
  return step;
}

solidStateEnthalpy()

double eonc::neb::solidStateEnthalpy	(	const Matter &	image,
		const Matter &	reference,
		double	pressure )

Potential energy plus P : (h0^{-1} (h-h0)) * V0.

Pressure is hydrostatic, in eV/Angstrom^3. A zero pressure returns the potential energy.

Definition at line 191 of file SolidStateNEB.cpp.

                                           {
  const double energy = image.getPotentialEnergy();
  if (pressure == 0.0) {
    return energy;
  }
  const Matrix3d h0 = reference.getCell();
  const Matrix3d strain = h0.inverse() * (image.getCell() - h0);
  const double volume = std::abs(h0.determinant());
  return energy + pressure * strain.trace() * volume;
}

solidStateJacobian()

double eonc::neb::solidStateJacobian	(	double	meanVolume,
		int	nAtoms,
		double	weight )

J = (V/N)^{1/3} * N^{1/2} * weight, with V the mean endpoint volume.

doi:10.1063/1.3684549

Definition at line 40 of file SolidStateNEB.cpp.

                                                                        {
  if (!(meanVolume > 0.0) || nAtoms < 1 || !(weight > 0.0)) {
    throw std::invalid_argument(
        "solid-state Jacobian needs a positive volume, atom count, and weight");
  }
  const double n = static_cast<double>(nAtoms);
  return std::cbrt(meanVolume / n) * std::sqrt(n) * weight;
}

writePathCon()

eonc::io::IoStatus eonc::neb::writePathCon ( const std::vector< std::shared_ptr< Matter > > & path, const std::vector< std::shared_ptr< AtomMatrix > > & tangent, const std::vector< std::shared_ptr< EigenmodeStrategy > > & eigenmode_solvers, long numImages, bool estimateEigenvalues, std::string filename, std::optional< size_t > bandIndex, double referenceEnergy )

nodiscard

Write a NEB band as a multi-frame .con via readcon ConFrameBuilder::clone().

Definition at line 374 of file NEBSplineExtrema.cpp.

                                                           {
  auto frames =
      pathToConFrames(path, tangent, eigenmode_solvers, numImages,
                      estimateEigenvalues, bandIndex, referenceEnergy);
  if (frames.empty()) {
    return eonc::io::IoStatus::InvalidArgument;
  }
  return eonc::io::writeConFrames(std::move(filename), frames);
}

zeroTranslation()

void eonc::neb::zeroTranslation	(	AtomMatrix &	projectedForce,
		int	nFreeAtoms,
		int	nAtoms )

Zero net translational force for fully free systems.

Definition at line 105 of file NEBForceProjection.cpp.

                                                                             {
  if (nFreeAtoms == nAtoms) {
    for (int j = 0; j <= 2; j++) {
      double translationMag = projectedForce.col(j).sum();
      int natoms = projectedForce.col(j).size();
      projectedForce.col(j).array() -=
          translationMag / static_cast<double>(natoms);
    }
  }
}