eonc::neb Namespace Reference

Classes
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	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)
	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)
	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.
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)
	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)
	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)
	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.

Typedef Documentation

ProjectionStrategy

using eonc::neb::ProjectionStrategy

Initial value:

 
std::variant<PlainEB, NEB_Projection, DNEB_Projection>

Definition at line 49 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{};
  }
  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 65 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{};
}

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 &	fPerp )

Compute the doubly-nudged elastic band perpendicular spring force.

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;
 
    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 )

Compute the DNEB force component for a given image.

Returned separately so it can be added to the CI force when DNEB is active.

Definition at line 47 of file NEBProjection.cpp.

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

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 74 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());
      tangentEndpoint.normalize();
      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());
      tangentEndpoint.normalize();
      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
      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 * 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;
}

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;
}

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 )

nodiscard

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

Empty on invalid path size. Does not write to disk.

Definition at line 251 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);
 
  for (long i = 0; i <= numImages + 1; i++) {
    if (i > 0) {
      distTotal += path[i]->distanceTo(*path[i - 1]);
    }
    metas.push_back(neb_frame_metadata(path, tangent, eigenmode_solvers,
                                       numImages, estimateEigenvalues, i,
                                       distTotal, bandIndex));
  }
 
  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 )

Print NEB image data to log and optionally to file.

Definition at line 166 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");
  }
 
  tangentStart.normalize();
  tangentEnd.normalize();
 
  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 = path[0]->getPotentialEnergy();
 
  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);
      }
    }
  }
}

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 )

nodiscard

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

Definition at line 279 of file NEBSplineExtrema.cpp.

                                   {
  auto frames = pathToConFrames(path, tangent, eigenmode_solvers, numImages,
                                estimateEigenvalues, bandIndex);
  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 103 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);
    }
  }
}