eonc::tunneling Namespace Reference

Classes
struct	Instanton
struct	InstantonOptions
class	Profile
	Energy along the path, interpolated as a monotone cubic (Fritsch and Carlson), flat at both ends because the ends of a band are minima. More...
struct	RateInstanton
struct	RateInstantonOptions
struct	RingSpectrum
	Spectrum of a closed ring's Hessian without forming it. More...
struct	Splitting

Typedefs
using	BatchPotential
	V (eV) and dV/dq (eV / (amu^0.5 Angstrom)) at every point of q, all in one call so a potential can spread the beads over its calculators.
using	BeadHessian = std::function<MatrixXd(long j, const VectorXd &q)>
	The mass-weighted Hessian d2V/dq2 at interior bead j (1..P-1).
using	RingBeadHessian = std::function<MatrixXd(long j, const VectorXd &q)>
	The mass-weighted Hessian d2V/dq2 at ring bead j (0..N-1).

Functions
double	massWeightedDistance (const Matter &a, const Matter &b)
	sqrt(sum_i m_i |b_i - a_i|^2) under the minimum image of a's cell.
std::vector< double >	massWeightedPath (const std::vector< std::shared_ptr< Matter > > &band)
	Cumulative mass-weighted arc length at each image of a band.
double	wellCurvature (const Profile &p, bool leftEnd)
	d2V/ds2 at one end of the path, in eV / (amu Angstrom^2), from a least squares fit of a s^2 + b s^3 to the points within half the barrier above that end.
double	hbarOmega (double curvature)
	hbar omega in eV for a mass-weighted curvature.
double	wkbAction (const Profile &p, double energy, int points=4001)
	(1/hbar) integral sqrt(2 (V(s) - E)) ds over the path where V > E.
Splitting	wkbSplitting (const Profile &p, double hwReactant, double hwProduct)
	delta0 = (hbar omega / pi) exp(-S) with the Landau and Lifshitz prefactor.
std::vector< VectorXd >	ringFromPath (const std::vector< VectorXd > &path, const std::vector< double > &energies, double betaHbar, long beads)
	Closed ring of beads samples of path whose imaginary-time period is betaHbar.
double	wkbLogRateAlongPath (const Profile &profile, double beta, double hwReactant)
	ln(k), k in 1/time, for the one-dimensional thermal rate along profile.
Splitting	bandSplitting (const std::vector< std::shared_ptr< Matter > > &band, double referenceEnergy)
	The splitting of a converged band, with the well frequencies from the band's curvature at each end.
double	pathOmega (const MatrixXd &hessStart, const MatrixXd &hessEnd, const VectorXd &start, const VectorXd &end)
	omega along the straight line between the minima from the curvature of each well there, the larger of the two; in 1 / time.
Instanton	optimizeInstanton (const VectorXd &start, const VectorXd &end, double betaHbar, std::vector< VectorXd > guess, const BatchPotential &potential, const InstantonOptions &options)
	Minimises the action from guess (P + 1 beads, ends at the minima, or empty for a tanh kink along the straight line).
void	instantonSplitting (Instanton &inst, const BeadHessian &hessian, const MatrixXd &hessStart, const MatrixXd &hessEnd)
	Fills delta0, zeroMode and modeSeparation from the bead Hessians and the Hessians of the two minima.
double	crossoverTemperature (const MatrixXd &hessSaddle)
	T_c = hbar omega_b / (2 pi kB) from the mass-weighted Hessian at the saddle, in K; throws when the Hessian has no negative eigenvalue.
RingSpectrum	ringSpectrum (const std::vector< MatrixXd > &beadHessians, double c, const std::vector< VectorXd > &tau)
	The ring Hessian of bead Hessians beadHessians (d2V/dq2 at each of the N beads) and spring constant c, with the normalised direction tau (N beads) projected out through the determinant lemma det(J + tau tau^T) = det' J when J tau = 0.
double	cyclicRingLogAbsDet (double c, const std::vector< MatrixXd > &diag)
	log|det| of the cyclic block-tridiagonal ring Hessian.
std::vector< VectorXd >	cyclicRingSolve (double c, const std::vector< MatrixXd > &diag, const std::vector< VectorXd > &rhs)
	Solves that same cyclic ring Hessian.
RateInstanton	optimizeRateInstanton (const VectorXd &saddle, const MatrixXd &hessSaddle, double beta, std::vector< VectorXd > guess, const BatchPotential &potential, const RateInstantonOptions &options)
	Finds the rate instanton at inverse temperature beta (1 / eV).
void	instantonRate (RateInstanton &inst, const RingBeadHessian &hessian, const MatrixXd &hessReactant, double vReactant, const MatrixXd &hessSaddle=MatrixXd(), double vSaddle=0.0, long rigidModes=0, long denseLimit=4096)
	Fills the rate from the bead Hessians, the reactant minimum's Hessian and energy, and optionally the saddle's Hessian and energy for the classical comparison (pass an empty matrix to skip it).
double	parabolicFactor (double temperature, double crossover)
	(pi T_c / T) / sin(pi T_c / T).
double	harmonicTstLogRate (const MatrixXd &hessReactant, const MatrixXd &hessSaddle, double beta, double barrier, long rigidModes)
	ln(k) for classical harmonic transition-state theory, k in 1/time.
double	quantumHarmonicTstLogRate (const MatrixXd &hessReactant, const MatrixXd &hessSaddle, double beta, double barrier, long rigidModes)
	ln(k) for quantum harmonic transition-state theory, k in 1/time: (1 / (2 pi beta hbar)) prod_r 2 sinh(beta hbar omega_r / 2) / prod'_s 2 sinh(beta hbar omega_s / 2) exp(-beta barrier), the zero-point and quantised partition functions of every bound mode; the saddle's unstable mode leaves the product.

Variables
constexpr double	kHbar = 0.06465415129579072
	hbar in eV^0.5 amu^0.5 Angstrom, correctly rounded from the exact SI values (h = 6.62607015e-34 J s, e = 1.602176634e-19 C) and the CODATA 2022 dalton, 1.66053906892e-27 kg.
constexpr double	kBoltzmann = 8.617333262145177e-5
	Boltzmann constant in eV / K, correctly rounded from the exact 1.380649e-23 J / K.
constexpr double	kTimeUnitSeconds = 1.0180505717871193e-14
	One unit of time, sqrt(amu Angstrom^2 / eV), in seconds.

Typedef Documentation

BatchPotential

using eonc::tunneling::BatchPotential

Initial value:

 
std::function<void(const std::vector<VectorXd> &q, std::vector<double> &v,
                   std::vector<VectorXd> &grad)>

V (eV) and dV/dq (eV / (amu^0.5 Angstrom)) at every point of q, all in one call so a potential can spread the beads over its calculators.

Definition at line 152 of file Tunneling.h.

BeadHessian

using eonc::tunneling::BeadHessian = std::function<MatrixXd(long j, const VectorXd &q)>

The mass-weighted Hessian d2V/dq2 at interior bead j (1..P-1).

Definition at line 157 of file Tunneling.h.

RingBeadHessian

using eonc::tunneling::RingBeadHessian = std::function<MatrixXd(long j, const VectorXd &q)>

The mass-weighted Hessian d2V/dq2 at ring bead j (0..N-1).

Definition at line 359 of file Tunneling.h.

Function Documentation

bandSplitting()

Splitting eonc::tunneling::bandSplitting	(	const std::vector< std::shared_ptr< Matter > > &	band,
		double	referenceEnergy )

The splitting of a converged band, with the well frequencies from the band's curvature at each end.

Definition at line 458 of file Tunneling.cpp.

                                                {
  std::vector<double> v;
  v.reserve(band.size());
  for (const auto &image : band) {
    v.push_back(image->getPotentialEnergy() - referenceEnergy);
  }
  const Profile p(massWeightedPath(band), std::move(v));
  return wkbSplitting(p, hbarOmega(wellCurvature(p, true)),
                      hbarOmega(wellCurvature(p, false)));
}

crossoverTemperature()

double eonc::tunneling::crossoverTemperature

(

const MatrixXd &

hessSaddle

)

T_c = hbar omega_b / (2 pi kB) from the mass-weighted Hessian at the saddle, in K; throws when the Hessian has no negative eigenvalue.

Definition at line 860 of file Tunneling.cpp.

                                                        {
  const Eigen::SelfAdjointEigenSolver<MatrixXd> es(
      0.5 * (hessSaddle + hessSaddle.transpose()));
  const double lambda = es.eigenvalues()(0);
  if (!(lambda < 0.0)) {
    throw std::invalid_argument(
        "crossoverTemperature: the saddle Hessian has no negative eigenvalue");
  }
  return kHbar * std::sqrt(-lambda) / (2.0 * std::numbers::pi * kBoltzmann);
}

cyclicRingLogAbsDet()

double eonc::tunneling::cyclicRingLogAbsDet	(	double	c,
		const std::vector< MatrixXd > &	diag )

log|det| of the cyclic block-tridiagonal ring Hessian.

Each diag[j] already contains the bead Hessian plus 2 c I, and the neighbour coupling is -c I, including the corner that closes the ring. A singular ring returns -infinity. Throws when the open chain is singular.

Definition at line 1580 of file Tunneling.cpp.

                                                                      {
  requireCyclicBlocks(c, diag, "cyclicRingLogAbsDet");
  return CyclicFactor(c, diag).logAbs;
}

cyclicRingSolve()

std::vector< VectorXd > eonc::tunneling::cyclicRingSolve	(	double	c,
		const std::vector< MatrixXd > &	diag,
		const std::vector< VectorXd > &	rhs )

Solves that same cyclic ring Hessian.

Throws when the ring is singular or the right-hand side does not match the blocks.

Definition at line 1585 of file Tunneling.cpp.

                                                                      {
  requireCyclicBlocks(c, diag, "cyclicRingSolve");
  if (static_cast<long>(rhs.size()) != static_cast<long>(diag.size())) {
    throw std::invalid_argument(
        "cyclicRingSolve: one right-hand side per bead");
  }
  for (const auto &row : rhs) {
    if (row.size() != diag.front().rows()) {
      throw std::invalid_argument(
          "cyclicRingSolve: the right-hand side does not match the blocks");
    }
  }
  return CyclicFactor(c, diag).solve(rhs);
}

harmonicTstLogRate()

double eonc::tunneling::harmonicTstLogRate	(	const MatrixXd &	hessReactant,
		const MatrixXd &	hessSaddle,
		double	beta,
		double	barrier,
		long	rigidModes )

ln(k) for classical harmonic transition-state theory, k in 1/time.

rigidModes eigenvalues nearest zero are omitted at each Hessian. The saddle's most negative eigenvalue is the barrier mode and leaves the product.

Definition at line 3151 of file Tunneling.cpp.

                                                           {
  if (!(beta > 0.0) || hessReactant.size() == 0 || hessSaddle.size() == 0 ||
      hessReactant.rows() != hessReactant.cols() ||
      hessSaddle.rows() != hessSaddle.cols() ||
      hessReactant.rows() != hessSaddle.rows() || rigidModes < 0) {
    throw std::invalid_argument(
        "harmonicTstLogRate: need beta > 0, matching square Hessians and a "
        "non-negative rigid-mode count");
  }
  const Eigen::SelfAdjointEigenSolver<MatrixXd> er(
      0.5 * (hessReactant + hessReactant.transpose()), Eigen::EigenvaluesOnly);
  const Eigen::SelfAdjointEigenSolver<MatrixXd> es(
      0.5 * (hessSaddle + hessSaddle.transpose()), Eigen::EigenvaluesOnly);
  const VectorXd &lr = er.eigenvalues();
  const VectorXd &ls = es.eigenvalues();
  if (!(ls(0) < 0.0)) {
    throw std::invalid_argument(
        "harmonicTstLogRate: the saddle Hessian has no negative eigenvalue");
  }
  const std::vector<bool> rigidR = nearestZero(lr, rigidModes);
  const std::vector<bool> rigidS = nearestZero(ls, rigidModes);
  double logRatio = 0.0;
  for (long m = 0; m < lr.size(); ++m) {
    if (!rigidR[static_cast<size_t>(m)]) {
      logRatio += 0.5 * std::log(lr(m));
    }
  }
  // Eigenvalue 0 is the unstable mode, the most negative.
  for (long m = 1; m < ls.size(); ++m) {
    if (!rigidS[static_cast<size_t>(m)]) {
      logRatio -= 0.5 * std::log(std::abs(ls(m)));
    }
  }
  return logRatio - std::log(2.0 * std::numbers::pi) - beta * barrier;
}

hbarOmega()

double eonc::tunneling::hbarOmega

(

double

curvature

)

hbar omega in eV for a mass-weighted curvature.

Definition at line 144 of file Tunneling.cpp.

                                   {
  if (!(curvature > 0.0)) {
    throw std::invalid_argument("a well needs a positive curvature");
  }
  return kHbar * std::sqrt(curvature);
}

instantonRate()

void eonc::tunneling::instantonRate	(	RateInstanton &	inst,
		const RingBeadHessian &	hessian,
		const MatrixXd &	hessReactant,
		double	vReactant,
		const MatrixXd &	hessSaddle = MatrixXd(),
		double	vSaddle = 0.0,
		long	rigidModes = 0,
		long	denseLimit = 4096 )

Fills the rate from the bead Hessians, the reactant minimum's Hessian and energy, and optionally the saddle's Hessian and energy for the classical comparison (pass an empty matrix to skip it).

rigidModes is the count of rigid-body zero modes to omit: the translations, plus a rotation only when the reactant Hessian leaves it null (a free cluster has them, a crystal does not, and an atom held fixed has none). They leave the centroid factors, so the rotational and translational partition functions of reactant and instanton cancel. Up to denseLimit ring degrees of freedom the product is the dense eigenproduct, checked against the cyclic block determinant. Beyond that, and for a limit of 0, the block determinant is used and the eigenvalues nearest zero come from inverse iteration on that factorisation. A negative limit forces the dense product.

Definition at line 2963 of file Tunneling.cpp.

                                    {
  const long N = static_cast<long>(inst.beads.size());
  if (N < 4 || !(inst.betaN > 0.0)) {
    throw std::invalid_argument("instantonRate: no optimised ring");
  }
  const long f = inst.beads.front().size();
  if (rigidModes < 0 || rigidModes > f) {
    throw std::invalid_argument(
        "instantonRate: rigidModes exceeds the degrees of freedom");
  }
  const double bnh = inst.betaN * kHbar;
  const double c = 1.0 / (bnh * bnh);
  const MatrixXd eye = MatrixXd::Identity(f, f);
 
  std::vector<MatrixXd> hBead(static_cast<size_t>(N));
  std::vector<MatrixXd> diag(static_cast<size_t>(N));
  for (long j = 0; j < N; ++j) {
    const MatrixXd h = hessian(j, inst.beads[static_cast<size_t>(j)]);
    if (h.rows() != f || h.cols() != f) {
      throw std::runtime_error("instantonRate: bead Hessian size");
    }
    hBead[static_cast<size_t>(j)] = 0.5 * (h + h.transpose());
    diag[static_cast<size_t>(j)] =
        hBead[static_cast<size_t>(j)] + 2.0 * c * eye;
  }
 
  const Eigen::SelfAdjointEigenSolver<MatrixXd> er(
      0.5 * (hessReactant + hessReactant.transpose()));
  const VectorXd &lr = er.eigenvalues();
  const std::vector<bool> rigidR = nearestZero(lr, rigidModes);
  MatrixXd nullBasis(f, 0);
  // The rigid vectors leave the product whether or not every bead Hessian
  // annihilates them exactly; a finite-difference Hessian never does, and
  // the reactant side drops the same modes at k = 0.
  for (long m = 0; m < lr.size(); ++m) {
    if (!rigidR[static_cast<size_t>(m)]) {
      continue;
    }
    nullBasis.conservativeResize(f, nullBasis.cols() + 1);
    nullBasis.col(nullBasis.cols() - 1) = er.eigenvectors().col(m);
  }
  if (rigidModes >= f) {
    throw std::runtime_error("instantonRate: every direction is a rigid mode");
  }
  // det' through the block chain: the cyclic zero mode and the rigid null
  // vectors leave the product by the determinant lemma, the inertia comes
  // from the Schur complements, and the N f by N f matrix is never formed.
  std::vector<VectorXd> cycle(static_cast<size_t>(N));
  double cycleNorm = 0.0;
  for (long j = 0; j < N; ++j) {
    cycle[static_cast<size_t>(j)] =
        0.5 * (inst.beads[static_cast<size_t>((j + 1) % N)] -
               inst.beads[static_cast<size_t>((j + N - 1) % N)]);
    cycleNorm += cycle[static_cast<size_t>(j)].squaredNorm();
  }
  if (!(cycleNorm > 0.0)) {
    throw std::runtime_error(
        "instantonRate: the beads coincide, so the ring has collapsed");
  }
  scale(cycle, 1.0 / std::sqrt(cycleNorm));
  // Each omitted direction is lifted by a spring-sized curvature c, far
  // above any physical near-zero eigenvalue, and the lift comes off the
  // log-determinant again: det(J + c u u^T) = c det' J when J u = 0.
  std::vector<std::vector<VectorXd>> dropped{cycle};
  std::vector<double> kappas{c};
  for (long r = 0; r < nullBasis.cols(); ++r) {
    dropped.emplace_back(
        static_cast<size_t>(N),
        (nullBasis.col(r) / std::sqrt(static_cast<double>(N))).eval());
    kappas.push_back(c);
  }
  const WoodburyRing ring(c, diag, true, dropped, kappas, true);
  if (!ring.ok() || !std::isfinite(ring.logAbsDet())) {
    throw std::runtime_error(
        "instantonRate: the ring Hessian is singular and the zero mode was "
        "not removed with the rigid modes");
  }
  inst.zeroEigenvalue = dot(cycle, applyDiagonal(diag, c, true, cycle));
  inst.negativeModes = ring.negative();
  // The lowest ring eigenvalue, for the report, from products alone.
  {
    auto applyFull = [&](const std::vector<VectorXd> &vec) {
      return applyDiagonal(diag, c, true, vec);
    };
    const long dim = N * f;
    const long steps = std::min(dim, static_cast<long>(60));
    std::vector<VectorXd> start = cycle;
    std::uint64_t h = 0x9E3779B97F4A7C15ULL;
    for (auto &bead : start) {
      for (long a2 = 0; a2 < bead.size(); ++a2) {
        h ^= h << 13;
        h ^= h >> 7;
        h ^= h << 17;
        bead(a2) += 0.1 * (static_cast<double>(h >> 11) * 0x1.0p-53 - 0.5);
      }
    }
    const std::vector<RingMode> modes =
        lowestRingModes(applyFull, std::move(start), steps);
    inst.negativeEigenvalue = 0.0;
    for (const auto &mode : modes) {
      if (mode.theta < inst.negativeEigenvalue &&
          std::abs(dot(mode.vector, cycle)) < 0.5) {
        inst.negativeEigenvalue = mode.theta;
      }
    }
    // A numerical null eigenvalue can sit just below zero off the cycle,
    // where the lift along tau does not reach it; it is not a second
    // unstable mode when it is tiny next to the barrier curvature.
    if (inst.negativeModes > 1 && inst.negativeEigenvalue < 0.0) {
      long tiny = 0;
      for (const auto &mode : modes) {
        if (mode.theta < 0.0 && mode.theta > 1e-3 * inst.negativeEigenvalue &&
            std::abs(dot(mode.vector, cycle)) < 0.5) {
          ++tiny;
        }
      }
      inst.negativeModes = std::max(1L, inst.negativeModes - tiny);
    }
  }
  const long nDrop = 1 + nullBasis.cols();
  const double logDetPrime =
      ring.logAbsDet() - static_cast<double>(nDrop) * std::log(c);
  const double logProd =
      static_cast<double>(N * f - nDrop) * std::log(bnh) + 0.5 * logDetPrime;
 
  inst.logRateTimesZr = -std::log(bnh) +
                        0.5 * std::log(inst.bN / (2.0 * std::numbers::pi *
                                                  inst.betaN * kHbar * kHbar)) -
                        logProd - inst.betaN * inst.ringPotential;
 
  for (long m = 0; m < lr.size(); ++m) {
    if (!rigidR[static_cast<size_t>(m)] && !(lr(m) > 0.0)) {
      throw std::runtime_error(
          "instantonRate: the reactant Hessian is not positive definite");
    }
  }
  // The rigid modes leave the centroid (k = 0) factor, as the ring's own
  // rigid modes leave its product; both carry them as free particles for
  // k > 0.
  double logZr = -inst.beta * vReactant;
  for (long k = 0; k < N; ++k) {
    const double sk = std::sin(std::numbers::pi * static_cast<double>(k) /
                               static_cast<double>(N));
    for (long m = 0; m < lr.size(); ++m) {
      if (k == 0 && rigidR[static_cast<size_t>(m)]) {
        continue;
      }
      const double l = rigidR[static_cast<size_t>(m)] ? 0.0 : lr(m);
      logZr -= std::log(bnh) + 0.5 * std::log(l + 4.0 * c * sk * sk);
    }
  }
  inst.logZr = logZr;
  inst.logRate = inst.logRateTimesZr - logZr;
  inst.rate = std::exp(inst.logRate) / kTimeUnitSeconds;
  inst.effectiveBarrier =
      -std::log(2.0 * std::numbers::pi * kHbar * inst.beta) / inst.beta -
      inst.logRate / inst.beta;
 
  if (hessSaddle.size() > 0) {
    inst.classicalLogRate = harmonicTstLogRate(
        hessReactant, hessSaddle, inst.beta, vSaddle - vReactant, rigidModes);
    inst.classicalRate = std::exp(inst.classicalLogRate) / kTimeUnitSeconds;
  }
}

instantonSplitting()

void eonc::tunneling::instantonSplitting	(	Instanton &	inst,
		const BeadHessian &	hessian,
		const MatrixXd &	hessStart,
		const MatrixXd &	hessEnd )

Fills delta0, zeroMode and modeSeparation from the bead Hessians and the Hessians of the two minima.

Definition at line 781 of file Tunneling.cpp.

                                                                            {
  const long P = static_cast<long>(inst.path.size()) - 1;
  if (P < 4 || !(inst.dtau > 0.0)) {
    throw std::invalid_argument("instantonSplitting: no optimised path");
  }
  const double dtau = inst.dtau;
  const double c = 1.0 / dtau;
  const long n = inst.path.front().size();
  const MatrixXd spring = 2.0 * c * MatrixXd::Identity(n, n);
 
  std::vector<MatrixXd> diag;
  diag.reserve(static_cast<size_t>(P - 1));
  for (long j = 1; j < P; ++j) {
    const MatrixXd h = hessian(j, inst.path[static_cast<size_t>(j)]);
    if (h.rows() != n || h.cols() != n) {
      throw std::runtime_error("instantonSplitting: bead Hessian size");
    }
    diag.push_back(spring + dtau * 0.5 * (h + h.transpose()));
  }
  const BlockChain chain(c, diag);
 
  auto wellLogDet = [&](const MatrixXd &h) {
    const std::vector<MatrixXd> d(static_cast<size_t>(P - 1),
                                  spring + dtau * 0.5 * (h + h.transpose()));
    const BlockChain well(c, d);
    if (well.sign() < 0) {
      throw std::runtime_error(
          "instantonSplitting: a well Hessian is not positive definite");
    }
    return well.logAbsDet();
  };
  const double logDetWell = 0.5 * (wellLogDet(hessStart) + wellLogDet(hessEnd));
 
  // The zero mode is the kink's translation in imaginary time, along the
  // discrete velocity v; det' J = det J (v^T J^-1 v) for v its eigenvector.
  std::vector<VectorXd> v(static_cast<size_t>(P - 1));
  for (long j = 1; j < P; ++j) {
    v[static_cast<size_t>(j - 1)] = inst.path[static_cast<size_t>(j + 1)] -
                                    inst.path[static_cast<size_t>(j - 1)];
  }
  scale(v, 1.0 / std::sqrt(dot(v, v)));
  const double vJv = dot(v, chain.solve(v));
  const int signPrime = chain.sign() * (vJv < 0.0 ? -1 : 1);
  if (signPrime < 0) {
    throw std::runtime_error(
        "instantonSplitting: the path is not a minimum of the action "
        "(a negative mode besides the kink's translation)");
  }
  inst.zeroMode = 1.0 / vJv;
  const double logDetPrime = chain.logAbsDet() + std::log(std::abs(vJv));
 
  // Next eigenvalue: inverse iteration orthogonal to v.
  std::vector<VectorXd> w(v.size());
  for (size_t k = 0; k < w.size(); ++k) {
    w[k].resize(n);
    for (long i = 0; i < n; ++i) {
      w[k](i) = std::sin(0.7 * static_cast<double>(k) +
                         1.3 * static_cast<double>(i) + 0.1);
    }
  }
  double lambda1 = 0.0;
  for (int it = 0; it < 40; ++it) {
    const double proj = dot(v, w);
    for (size_t k = 0; k < w.size(); ++k) {
      w[k] -= proj * v[k];
    }
    scale(w, 1.0 / std::sqrt(dot(w, w)));
    std::vector<VectorXd> z = chain.solve(w);
    lambda1 = 1.0 / dot(w, z);
    w = std::move(z);
  }
  inst.modeSeparation = std::abs(lambda1 / inst.zeroMode);
 
  inst.delta0 = 2.0 * kHbar *
                std::sqrt(inst.s0 / (2.0 * std::numbers::pi * kHbar * dtau)) *
                std::exp(0.5 * (logDetWell - logDetPrime) - inst.action);
}

massWeightedDistance()

double eonc::tunneling::massWeightedDistance	(	const Matter &	a,
		const Matter &	b )

sqrt(sum_i m_i |b_i - a_i|^2) under the minimum image of a's cell.

Definition at line 34 of file Tunneling.cpp.

                                                              {
  if (a.numberOfAtoms() != b.numberOfAtoms()) {
    throw std::invalid_argument("the structures hold different atom counts");
  }
  const AtomMatrix dr = a.pbc(b.getPositions() - a.getPositions());
  const auto mass = a.getMasses();
  double sum = 0.0;
  for (long i = 0; i < a.numberOfAtoms(); ++i) {
    if (mass(i) <= 0.0) {
      throw std::invalid_argument(
          "every atom needs a positive mass for a mass-weighted path");
    }
    sum += mass(i) * dr.row(i).squaredNorm();
  }
  return std::sqrt(sum);
}

massWeightedPath()

std::vector< double > eonc::tunneling::massWeightedPath

(

const std::vector< std::shared_ptr< Matter > > &

band

)

Cumulative mass-weighted arc length at each image of a band.

Definition at line 52 of file Tunneling.cpp.

                                                               {
  std::vector<double> s{0.0};
  s.reserve(band.size());
  for (size_t i = 1; i < band.size(); ++i) {
    s.push_back(s.back() + massWeightedDistance(*band[i - 1], *band[i]));
  }
  return s;
}

optimizeInstanton()

Instanton eonc::tunneling::optimizeInstanton	(	const VectorXd &	start,
		const VectorXd &	end,
		double	betaHbar,
		std::vector< VectorXd >	guess,
		const BatchPotential &	potential,
		const InstantonOptions &	options )

Minimises the action from guess (P + 1 beads, ends at the minima, or empty for a tanh kink along the straight line).

Evaluates the interior beads once per iteration and line-search step.

Definition at line 633 of file Tunneling.cpp.

                                                             {
  const long P = options.beads;
  if (P < 4 || !(betaHbar > 0.0) || start.size() != end.size()) {
    throw std::invalid_argument("optimizeInstanton: need P >= 4, beta hbar > 0 "
                                "and ends of one dimension");
  }
  Instanton inst;
  inst.betaHbar = betaHbar;
  inst.dtau = betaHbar / static_cast<double>(P);
  const double dtau = inst.dtau;
 
  std::vector<double> vEnds;
  std::vector<VectorXd> gEnds;
  potential({start, end}, vEnds, gEnds);
  if (vEnds.size() != 2) {
    throw std::runtime_error("instanton: potential returned the wrong count");
  }
  inst.asymmetry = vEnds[1] - vEnds[0];
 
  // Beads along the guess (or the straight line) on a tanh kink centred at
  // beta hbar / 2 whose width follows the harmonic decay of a well.
  if (guess.size() < 2) {
    guess = {start, end};
  }
  std::vector<double> cum(guess.size(), 0.0);
  for (size_t k = 1; k < guess.size(); ++k) {
    cum[k] = cum[k - 1] + (guess[k] - guess[k - 1]).norm();
  }
  if (!(cum.back() > 0.0)) {
    throw std::invalid_argument("optimizeInstanton: the two minima coincide");
  }
  const double width = betaHbar / (2.0 * options.betaHbarOmega);
  std::vector<VectorXd> x(static_cast<size_t>(P - 1));
  for (long j = 1; j < P; ++j) {
    const double tau = static_cast<double>(j) * dtau - 0.5 * betaHbar;
    const double f = 0.5 * (1.0 + std::tanh(tau / width));
    x[static_cast<size_t>(j - 1)] = alongPolyline(guess, cum, f);
  }
 
  // L-BFGS with a backtracking Armijo line search.
  ActionEval cur =
      evaluateAction(x, start, end, vEnds[0], vEnds[1], dtau, potential);
  std::deque<std::pair<std::vector<VectorXd>, std::vector<VectorXd>>> pairs;
  for (long it = 0; it < options.maxIterations; ++it) {
    inst.iterations = it;
    if (largestBeadNorm(cur.grad) / dtau < options.forceTolerance) {
      inst.converged = true;
      break;
    }
    std::vector<VectorXd> q = cur.grad;
    std::vector<double> alpha(pairs.size());
    for (size_t i = pairs.size(); i-- > 0;) {
      const double rho = 1.0 / dot(pairs[i].second, pairs[i].first);
      alpha[i] = rho * dot(pairs[i].first, q);
      for (size_t k = 0; k < q.size(); ++k) {
        q[k] -= alpha[i] * pairs[i].second[k];
      }
    }
    // Without history, half the inverse spring stiffness 2 / dtau.
    double gamma = dtau / 4.0;
    if (!pairs.empty()) {
      gamma = dot(pairs.back().first, pairs.back().second) /
              dot(pairs.back().second, pairs.back().second);
    }
    scale(q, gamma);
    for (size_t i = 0; i < pairs.size(); ++i) {
      const double rho = 1.0 / dot(pairs[i].second, pairs[i].first);
      const double beta = rho * dot(pairs[i].second, q);
      for (size_t k = 0; k < q.size(); ++k) {
        q[k] += (alpha[i] - beta) * pairs[i].first[k];
      }
    }
    // q is now the inverse-Hessian estimate times the gradient; step -q.
    double slope = -dot(cur.grad, q);
    if (!(slope < 0.0)) {
      pairs.clear();
      q = cur.grad;
      scale(q, dtau / 4.0);
      slope = -dot(cur.grad, q);
    }
    double step = 1.0;
    ActionEval next;
    std::vector<VectorXd> trial(x.size());
    bool accepted = false;
    for (int ls = 0; ls < 30; ++ls) {
      for (size_t k = 0; k < x.size(); ++k) {
        trial[k] = x[k] - step * q[k];
      }
      next = evaluateAction(trial, start, end, vEnds[0], vEnds[1], dtau,
                            potential);
      // Near the minimum the action changes by less than its round-off;
      // there a step that shrinks the gradient is progress too.
      const bool armijo = next.action <= cur.action + 1e-4 * step * slope;
      const bool flat = std::abs(next.action - cur.action) <=
                        1e-13 * std::max(1.0, std::abs(cur.action));
      if (armijo ||
          (flat && largestBeadNorm(next.grad) < largestBeadNorm(cur.grad))) {
        accepted = true;
        break;
      }
      step *= 0.5;
    }
    if (!accepted) {
      break;
    }
    std::vector<VectorXd> sk(x.size()), yk(x.size());
    for (size_t k = 0; k < x.size(); ++k) {
      sk[k] = trial[k] - x[k];
      yk[k] = next.grad[k] - cur.grad[k];
    }
    if (dot(sk, yk) > 0.0) {
      pairs.emplace_back(std::move(sk), std::move(yk));
      if (static_cast<long>(pairs.size()) > options.memory) {
        pairs.pop_front();
      }
    }
    x = std::move(trial);
    cur = std::move(next);
  }
  if (!inst.converged &&
      largestBeadNorm(cur.grad) / dtau < options.forceTolerance) {
    inst.converged = true;
  }
 
  inst.path.reserve(static_cast<size_t>(P + 1));
  inst.path.push_back(start);
  inst.path.insert(inst.path.end(), x.begin(), x.end());
  inst.path.push_back(end);
  inst.energies.reserve(static_cast<size_t>(P + 1));
  inst.energies.push_back(vEnds[0]);
  inst.energies.insert(inst.energies.end(), cur.v.begin(), cur.v.end());
  inst.energies.push_back(vEnds[1]);
  const double sWell = betaHbar * 0.5 * (vEnds[0] + vEnds[1]);
  inst.action = (cur.action - sWell) / kHbar;
  double s0 = 0.0;
  for (long j = 0; j < P; ++j) {
    s0 += (inst.path[static_cast<size_t>(j + 1)] -
           inst.path[static_cast<size_t>(j)])
              .squaredNorm();
  }
  inst.s0 = s0 / dtau;
  inst.symmetricEnough = std::abs(inst.asymmetry) * betaHbar / kHbar < 0.1;
  return inst;
}

optimizeRateInstanton()

RateInstanton eonc::tunneling::optimizeRateInstanton	(	const VectorXd &	saddle,
		const MatrixXd &	hessSaddle,
		double	beta,
		std::vector< VectorXd >	guess,
		const BatchPotential &	potential,
		const RateInstantonOptions &	options )

Finds the rate instanton at inverse temperature beta (1 / eV).

guess holds N beads, or is empty for a ring stretched along the saddle's unstable mode to where V has dropped by (1 - T / T_c) of the lower of the two barriers. At or below newtonLimit active coordinates the step is an index-1 Newton step on a Bofill Hessian, and below 0.75 T_c an empty guess cools from 0.85 T_c. Beyond that limit the step is minimum-mode following. An even ring whose beads match under j -> N - j is evaluated from one turning point to the other and mirrored. saddle and hessSaddle are mass-weighted.

Definition at line 2674 of file Tunneling.cpp.

                                                                         {
  const long N = options.beads;
  if (N < 4 || !(beta > 0.0) || hessSaddle.rows() != saddle.size()) {
    throw std::invalid_argument(
        "optimizeRateInstanton: need N >= 4, beta > 0 and a saddle Hessian "
        "of the saddle's dimension");
  }
  RateInstanton inst;
  inst.beta = beta;
  inst.betaN = beta / static_cast<double>(N);
  inst.temperature = 1.0 / (kBoltzmann * beta);
  inst.crossover = crossoverTemperature(hessSaddle);
  if (!(inst.temperature < inst.crossover)) {
    throw std::invalid_argument(
        "optimizeRateInstanton: T is at or above the crossover temperature; "
        "the ring collapses onto the saddle and steepest descent needs the "
        "parabolic barrier correction, of which classical transition-state "
        "theory is only the one-bead limit");
  }
  bool mirror = options.halfRing && N % 2 == 0;
  if (mirror && static_cast<long>(guess.size()) == N) {
    const long mid = N / 2;
    for (long j = 1; j < mid && mirror; ++j) {
      if ((guess[static_cast<size_t>(j)] - guess[static_cast<size_t>(N - j)])
              .norm() > 1e-8) {
        mirror = false;
      }
    }
  }
  const long active = mirror ? (N / 2 + 1) : N;
  if (options.newtonLimit > 0 &&
      active * saddle.size() <= options.newtonLimit) {
    return optimizeRateByNewton(saddle, hessSaddle, beta, std::move(guess),
                                potential, options);
  }
  const double bnh = inst.betaN * kHbar;
  const double c = 1.0 / (bnh * bnh);
 
  const ColMajorXd saddleCurvature =
      0.5 * (hessSaddle + hessSaddle.transpose());
  const Eigen::SelfAdjointEigenSolver<ColMajorXd> es(saddleCurvature);
  const VectorXd dir = es.eigenvectors().col(0);
 
  if (static_cast<long>(guess.size()) != N) {
    std::vector<double> v0;
    std::vector<VectorXd> g0;
    potential({saddle}, v0, g0);
    const double vS = v0.at(0);
    // Scan in steps of a quarter of the length over which the barrier's
    // curvature drops V by kB T_c.
    const double h = 0.25 * std::sqrt(2.0 * kBoltzmann * inst.crossover /
                                      -es.eigenvalues()(0));
    const long pts = 200;
    const double dPlus = sideDrop(saddle, dir, vS, 1.0, h, pts, potential);
    const double dMinus = sideDrop(saddle, dir, vS, -1.0, h, pts, potential);
    const double dMin = std::min(dPlus, dMinus);
    const double drop =
        (1.0 - inst.temperature / inst.crossover) *
        (std::isfinite(dMin) ? dMin : kBoltzmann * inst.crossover);
    const double sPlus =
        turningDistance(saddle, dir, vS, drop, 1.0, h, pts, potential);
    const double sMinus =
        turningDistance(saddle, dir, vS, drop, -1.0, h, pts, potential);
    guess.resize(static_cast<size_t>(N));
    for (long j = 0; j < N; ++j) {
      const double ct =
          std::cos(2.0 * std::numbers::pi * static_cast<double>(j) /
                   static_cast<double>(N));
      guess[static_cast<size_t>(j)] =
          saddle + dir * (ct >= 0.0 ? sPlus * ct : sMinus * ct);
    }
  }
 
  // An even count whose beads already match under j -> N-j is the
  // out-and-back instanton. The images are assigned equal, so the
  // potential on one half is copied, and the step stays on the closed ring.
  const bool wantMirror = options.halfRing && N % 2 == 0;
  std::vector<VectorXd> x = std::move(guess);
  bool fold = false;
  if (wantMirror) {
    fold = true;
    const long m = N / 2;
    for (long j = 1; j < m; ++j) {
      if ((x[static_cast<size_t>(j)] - x[static_cast<size_t>(N - j)]).norm() >
          1e-8) {
        fold = false;
        break;
      }
    }
  }
  BatchPotential evalPot = potential;
  if (fold) {
    evalPot = [&](const std::vector<VectorXd> &q, std::vector<double> &v,
                  std::vector<VectorXd> &g) {
      const long m = N / 2;
      bool sym = true;
      for (long j = 1; j < m; ++j) {
        if ((q[static_cast<size_t>(j)] - q[static_cast<size_t>(N - j)])
                .squaredNorm() != 0.0) {
          sym = false;
          break;
        }
      }
      if (!sym) {
        potential(q, v, g);
        return;
      }
      std::vector<VectorXd> uniq(static_cast<size_t>(m + 1));
      for (long j = 0; j <= m; ++j) {
        uniq[static_cast<size_t>(j)] = q[static_cast<size_t>(j)];
      }
      std::vector<double> vu;
      std::vector<VectorXd> gu;
      potential(uniq, vu, gu);
      v.assign(static_cast<size_t>(N), 0.0);
      g.assign(static_cast<size_t>(N), VectorXd());
      for (long j = 0; j <= m; ++j) {
        v[static_cast<size_t>(j)] = vu[static_cast<size_t>(j)];
        g[static_cast<size_t>(j)] = gu[static_cast<size_t>(j)];
      }
      for (long j = 1; j < m; ++j) {
        v[static_cast<size_t>(N - j)] = vu[static_cast<size_t>(j)];
        g[static_cast<size_t>(N - j)] = gu[static_cast<size_t>(j)];
      }
    };
  }
  auto symmetrize = [&](std::vector<VectorXd> &q) {
    if (!fold) {
      return;
    }
    const long m = N / 2;
    for (long j = 1; j < m; ++j) {
      const size_t a = static_cast<size_t>(j);
      const size_t b = static_cast<size_t>(N - j);
      const VectorXd mid = 0.5 * (q[a] + q[b]);
      q[a] = mid;
      q[b] = mid;
    }
  };
  symmetrize(x);
 
  RingEval cur = evaluateRing(x, c, evalPot, options.energyShift);
  // The unstable mode of the ring starts as every bead moving along the
  // saddle's unstable direction.
  std::vector<VectorXd> mode(x.size(), dir);
  double curvature = lowestMode(x, cur, c, evalPot, mode, options.lanczosFirst,
                                options.lanczosStep, fold);
  std::deque<std::pair<std::vector<VectorXd>, std::vector<VectorXd>>> pairs;
  auto effective = [&](const std::vector<VectorXd> &g) {
    const double par = dot(g, mode);
    std::vector<VectorXd> e = g;
    const double f = curvature < 0.0 ? 2.0 : 1.0;
    for (size_t j = 0; j < e.size(); ++j) {
      e[j] -= f * par * mode[j];
      if (!(curvature < 0.0)) {
        e[j] = -par * mode[j]; // climb along the mode only
      }
    }
    return e;
  };
  std::vector<VectorXd> geff = effective(cur.grad);
  long limit = options.maxIterations;
  for (long it = 0; it < limit; ++it) {
    inst.iterations = it;
    if (curvature < 0.0 && largestBeadNorm(cur.grad) < options.forceTolerance) {
      std::vector<VectorXd> oddMode;
      const double oddCurv =
          fold ? lowestOddMode(x, cur, c, potential, oddMode,
                               options.lanczosFirst, options.lanczosStep)
               : 0.0;
      if (!(oddCurv < -1e-3 * std::abs(curvature))) {
        inst.converged = true;
        break;
      }
      // A mirror-symmetric stationary point with a second unstable mode
      // odd under the mirror, such as two copies of the instanton on one
      // ring. The rest of the search runs on the whole ring from a kick
      // along that mode which lowers beta_N U_N by one in the quadratic
      // model.
      fold = false;
      evalPot = potential;
      const double kick = std::sqrt(2.0 / (inst.betaN * -oddCurv));
      for (size_t k = 0; k < x.size(); ++k) {
        x[k] += kick * oddMode[k];
      }
      cur = evaluateRing(x, c, evalPot, options.energyShift);
      curvature = lowestMode(x, cur, c, evalPot, mode, options.lanczosFirst,
                             options.lanczosStep, fold);
      pairs.clear();
      geff = effective(cur.grad);
      limit = it + options.maxIterations;
    }
    std::vector<VectorXd> trial(x.size());
    std::vector<VectorXd> d = geff;
    std::vector<double> alpha(pairs.size());
    for (size_t i = pairs.size(); i-- > 0;) {
      const double rho = 1.0 / dot(pairs[i].second, pairs[i].first);
      alpha[i] = rho * dot(pairs[i].first, d);
      for (size_t k = 0; k < d.size(); ++k) {
        d[k] -= alpha[i] * pairs[i].second[k];
      }
    }
    double gamma = 1.0 / (4.0 * c); // spring stiffness sets the first scale
    if (!pairs.empty()) {
      gamma = dot(pairs.back().first, pairs.back().second) /
              dot(pairs.back().second, pairs.back().second);
    }
    scale(d, gamma);
    for (size_t i = 0; i < pairs.size(); ++i) {
      const double rho = 1.0 / dot(pairs[i].second, pairs[i].first);
      const double b = rho * dot(pairs[i].second, d);
      for (size_t k = 0; k < d.size(); ++k) {
        d[k] += (alpha[i] - b) * pairs[i].first[k];
      }
    }
    if (!(dot(d, geff) > 0.0)) { // not a descent direction on geff
      pairs.clear();
      d = geff;
      scale(d, 1.0 / (4.0 * c));
    }
    const double big = largestBeadNorm(d);
    if (big > options.maxStep) {
      scale(d, options.maxStep / big);
    }
    for (size_t k = 0; k < x.size(); ++k) {
      trial[k] = x[k] - d[k];
    }
    symmetrize(trial);
    RingEval next = evaluateRing(trial, c, evalPot, options.energyShift);
    const double prevCurv = curvature;
    const long restart = options.lanczosRestart;
    curvature = lowestMode(trial, next, c, evalPot, mode, restart,
                           options.lanczosStep, fold);
    std::vector<VectorXd> geffNext = effective(next.grad);
    if ((prevCurv < 0.0) != (curvature < 0.0)) {
      pairs.clear();
    } else {
      std::vector<VectorXd> sk(x.size()), yk(x.size());
      for (size_t k = 0; k < x.size(); ++k) {
        sk[k] = trial[k] - x[k];
        yk[k] = geffNext[k] - geff[k];
      }
      if (dot(sk, yk) > 0.0) {
        pairs.emplace_back(std::move(sk), std::move(yk));
        if (static_cast<long>(pairs.size()) > options.memory) {
          pairs.pop_front();
        }
      }
    }
    x = std::move(trial);
    cur = std::move(next);
    geff = std::move(geffNext);
  }
  if (!inst.converged && curvature < 0.0 &&
      largestBeadNorm(cur.grad) < options.forceTolerance) {
    inst.converged = true;
  }
  inst.beads = x;
  inst.energies = cur.v;
  inst.ringPotential = cur.u;
  inst.bN = 0.0;
  for (size_t j = 0; j < x.size(); ++j) {
    inst.bN += (x[(j + 1) % x.size()] - x[j]).squaredNorm();
  }
  return inst;
}

parabolicFactor()

double eonc::tunneling::parabolicFactor	(	double	temperature,
		double	crossover )

(pi T_c / T) / sin(pi T_c / T).

Defined only above T_c. The factor diverges as T approaches T_c from above and tends to 1 at high T.

Definition at line 3131 of file Tunneling.cpp.

                                                             {
  if (!(temperature > 0.0) || !(crossover > 0.0)) {
    throw std::invalid_argument(
        "parabolicFactor: temperature and crossover must be positive");
  }
  if (!(temperature > crossover)) {
    throw std::invalid_argument(
        "parabolicFactor: T is at or below the crossover; the factor "
        "diverges there");
  }
  // beta hbar omega_b / 2 = pi T_c / T, since T_c = hbar omega_b / (2 pi kB).
  const double phase = std::numbers::pi * crossover / temperature;
  const double s = std::sin(phase);
  if (!(s > 0.0)) {
    throw std::invalid_argument(
        "parabolicFactor: the sine of the barrier phase is not positive");
  }
  return phase / s;
}

pathOmega()

double eonc::tunneling::pathOmega	(	const MatrixXd &	hessStart,
		const MatrixXd &	hessEnd,
		const VectorXd &	start,
		const VectorXd &	end )

omega along the straight line between the minima from the curvature of each well there, the larger of the two; in 1 / time.

Definition at line 622 of file Tunneling.cpp.

                                                             {
  const VectorXd d = (end - start).normalized();
  const double k = std::max(d.dot(hessStart * d), d.dot(hessEnd * d));
  if (!(k > 0.0)) {
    throw std::invalid_argument(
        "pathOmega: no positive curvature along the path at either minimum");
  }
  return std::sqrt(k);
}

quantumHarmonicTstLogRate()

double eonc::tunneling::quantumHarmonicTstLogRate	(	const MatrixXd &	hessReactant,
		const MatrixXd &	hessSaddle,
		double	beta,
		double	barrier,
		long	rigidModes )

ln(k) for quantum harmonic transition-state theory, k in 1/time: (1 / (2 pi beta hbar)) prod_r 2 sinh(beta hbar omega_r / 2) / prod'_s 2 sinh(beta hbar omega_s / 2) exp(-beta barrier), the zero-point and quantised partition functions of every bound mode; the saddle's unstable mode leaves the product.

Times the parabolic factor it is the N -> infinity ring-polymer rate above T_c, so it joins the instanton rate at the crossover; its high-temperature limit is harmonicTstLogRate.

Definition at line 3189 of file Tunneling.cpp.

                                                                  {
  if (!(beta > 0.0) || hessReactant.size() == 0 || hessSaddle.size() == 0 ||
      hessReactant.rows() != hessReactant.cols() ||
      hessSaddle.rows() != hessSaddle.cols() ||
      hessReactant.rows() != hessSaddle.rows() || rigidModes < 0) {
    throw std::invalid_argument(
        "quantumHarmonicTstLogRate: need beta > 0, matching square Hessians "
        "and a non-negative rigid-mode count");
  }
  const ColMajorXd hr = 0.5 * (hessReactant + hessReactant.transpose());
  const ColMajorXd hsd = 0.5 * (hessSaddle + hessSaddle.transpose());
  const Eigen::SelfAdjointEigenSolver<ColMajorXd> er(hr,
                                                     Eigen::EigenvaluesOnly);
  const Eigen::SelfAdjointEigenSolver<ColMajorXd> es(hsd,
                                                     Eigen::EigenvaluesOnly);
  const VectorXd lr = er.eigenvalues();
  const VectorXd ls = es.eigenvalues();
  if (!(ls(0) < 0.0)) {
    throw std::invalid_argument("quantumHarmonicTstLogRate: the saddle "
                                "Hessian has no negative eigenvalue");
  }
  const std::vector<bool> rigidR = nearestZero(lr, rigidModes);
  const std::vector<bool> rigidS = nearestZero(ls, rigidModes);
  const double bh = beta * kHbar;
  // ln(2 sinh(x / 2)) without overflow for large x.
  auto logTwoSinhHalf = [](double x) {
    return 0.5 * x + std::log1p(-std::exp(-x));
  };
  double logRatio = 0.0;
  for (long m = 0; m < lr.size(); ++m) {
    if (!rigidR[static_cast<size_t>(m)]) {
      logRatio += logTwoSinhHalf(bh * std::sqrt(lr(m)));
    }
  }
  for (long m = 1; m < ls.size(); ++m) {
    if (!rigidS[static_cast<size_t>(m)]) {
      logRatio -= logTwoSinhHalf(bh * std::sqrt(std::abs(ls(m))));
    }
  }
  return logRatio - std::log(2.0 * std::numbers::pi * bh) - beta * barrier;
}

ringFromPath()

std::vector< VectorXd > eonc::tunneling::ringFromPath	(	const std::vector< VectorXd > &	path,
		const std::vector< double > &	energies,
		double	betaHbar,
		long	beads )

Closed ring of beads samples of path whose imaginary-time period is betaHbar.

Bead 0 is the reactant-side turning point and bead N/2 the other; bead N - j repeats bead j. Throws when the path has no barrier, or when the period at the barrier top already exceeds betaHbar.

Definition at line 270 of file Tunneling.cpp.

                                                                {
  if (path.size() < 3 || path.size() != energies.size() || beads < 4 ||
      !(betaHbar > 0.0)) {
    throw std::invalid_argument(
        "ringFromPath: a path of at least three points with energies, "
        "N >= 4 and beta hbar > 0");
  }
  const long width = path.front().size();
  for (const auto &q : path) {
    if (q.size() != width) {
      throw std::invalid_argument("ringFromPath: the path changes dimension");
    }
  }
  const std::vector<double> s = arcLengths(path);
  const Profile profile(s, energies);
  double sTop = s.front();
  double vTop = -std::numeric_limits<double>::infinity();
  const int grid = 4000;
  for (int k = 0; k <= grid; ++k) {
    const double sk =
        s.front() + (s.back() - s.front()) * static_cast<double>(k) / grid;
    if (profile(sk) > vTop) {
      vTop = profile(sk);
      sTop = sk;
    }
  }
  const double vLow = std::max(energies.front(), energies.back());
  if (!(vTop > vLow)) {
    throw std::invalid_argument("ringFromPath: the path has no barrier");
  }
  auto period = [&](double energy) {
    const auto [sMinus, sPlus] = turningPoints(profile, sTop, energy);
    return 2.0 * halfPeriod(profile, sMinus, sPlus, energy);
  };
  // The crossover along the path from a parabola through the three input
  // points around the barrier top, not from the interpolant, whose slope is
  // clamped to zero at the top node.
  {
    size_t top = 0;
    for (size_t k = 1; k < energies.size(); ++k) {
      if (energies[k] > energies[top]) {
        top = k;
      }
    }
    if (top == 0 || top + 1 == energies.size()) {
      throw std::invalid_argument(
          "ringFromPath: the barrier top is an end of the path");
    }
    const double h1 = s[top] - s[top - 1], h2 = s[top + 1] - s[top];
    const double curvature =
        2.0 *
        (h1 * energies[top + 1] - (h1 + h2) * energies[top] +
         h2 * energies[top - 1]) /
        (h1 * h2 * (h1 + h2));
    if (!(curvature < 0.0)) {
      throw std::invalid_argument("ringFromPath: no curvature at the top");
    }
    const double tc = kHbar * std::sqrt(-curvature) / (2.0 * std::numbers::pi);
    if (!(kHbar / betaHbar < tc)) {
      throw std::invalid_argument(
          "ringFromPath: the temperature is at or above the crossover along "
          "this path");
    }
  }
  // Bracket the orbit energy geometrically above the lower end: the period
  // grows only logarithmically as E approaches a well bottom.
  const double span = vTop - vLow;
  double eHi = vTop - 1e-9 * span;
  double eLo = vLow + 1e-14 * span;
  if (period(eLo) < betaHbar) {
    // The path does not reach a long enough orbit. The lowest one it
    // holds is the start.
    eHi = eLo;
  }
  for (int k = 0; k < 200 && eHi > eLo; ++k) {
    const double e = vLow + std::sqrt((eLo - vLow) * (eHi - vLow));
    if (period(e) > betaHbar) {
      eLo = e;
    } else {
      eHi = e;
    }
    if (eHi - eLo < 1e-15 * span) {
      break;
    }
  }
  const double energy = 0.5 * (eLo + eHi);
  const auto [sMinus, sPlus] = turningPoints(profile, sTop, energy);
  std::vector<double> tau;
  std::vector<double> pos;
  const double half = halfPeriod(profile, sMinus, sPlus, energy, &tau, &pos);
  if (tau.size() < 2 || tau.size() != pos.size()) {
    throw std::invalid_argument("ringFromPath: the orbit has no length");
  }
  // Bead j sits at imaginary time j * beta hbar / N on the way from the
  // reactant-side turning point to the other side. The return repeats it.
  std::vector<VectorXd> ring(static_cast<size_t>(beads), VectorXd::Zero(width));
  for (long j = 0; j <= beads / 2; ++j) {
    const double t =
        std::min(half, half * 2.0 * static_cast<double>(j) / beads);
    const auto it = std::upper_bound(tau.begin(), tau.end(), t);
    const size_t k = std::clamp<size_t>(
        static_cast<size_t>(std::distance(tau.begin(), it)), 1, tau.size() - 1);
    const double seg = tau[k] - tau[k - 1];
    const double w = seg > 0.0 ? (t - tau[k - 1]) / seg : 0.0;
    const double sj =
        pos[k - 1] + std::clamp(w, 0.0, 1.0) * (pos[k] - pos[k - 1]);
    ring[static_cast<size_t>(j)] = atArcLength(path, s, sj);
    if (j > 0 && j < beads - j) {
      ring[static_cast<size_t>(beads - j)] = ring[static_cast<size_t>(j)];
    }
  }
  return ring;
}

ringSpectrum()

RingSpectrum eonc::tunneling::ringSpectrum	(	const std::vector< MatrixXd > &	beadHessians,
		double	c,
		const std::vector< VectorXd > &	tau )

The ring Hessian of bead Hessians beadHessians (d2V/dq2 at each of the N beads) and spring constant c, with the normalised direction tau (N beads) projected out through the determinant lemma det(J + tau tau^T) = det' J when J tau = 0.

Block LU of the open chain plus a low-rank correction for the closure and tau, O(N f^3).

Definition at line 1562 of file Tunneling.cpp.

                                                          {
  if (beadHessians.empty() || beadHessians.size() != tau.size()) {
    throw std::invalid_argument(
        "ringSpectrum: N bead Hessians and N tau blocks");
  }
  const std::vector<MatrixXd> diag = ringDiagonal(beadHessians, c, false);
  const WoodburyRing ring(c, diag, true, {tau}, {1.0}, true);
  if (!ring.ok()) {
    throw std::runtime_error("ringSpectrum: singular ring");
  }
  RingSpectrum out;
  out.logDetPrime = ring.logAbsDet();
  out.negativeModes = ring.negative();
  out.zeroEigenvalue = dot(tau, applyDiagonal(diag, c, true, tau));
  return out;
}

wellCurvature()

double eonc::tunneling::wellCurvature	(	const Profile &	p,
		bool	leftEnd )

d2V/ds2 at one end of the path, in eV / (amu Angstrom^2), from a least squares fit of a s^2 + b s^3 to the points within half the barrier above that end.

The cubic term takes the anharmonic part a parabola alone would fold into the curvature.

Definition at line 104 of file Tunneling.cpp.

                                                     {
  const auto &s = p.s();
  const auto &v = p.v();
  const size_t n = s.size();
  const double top = *std::max_element(v.begin(), v.end());
  const double floor = leftEnd ? v.front() : v.back();
  const double half = 0.5 * (top - floor);
  // Sums for the normal equations of y = a x^2 + b x^3.
  double s44 = 0, s45 = 0, s55 = 0, sy2 = 0, sy3 = 0;
  size_t used = 0;
  for (size_t j = 1; j < n; ++j) {
    const size_t i = leftEnd ? j : n - 1 - j;
    const double x = leftEnd ? s[i] - s.front() : s.back() - s[i];
    const double y = v[i] - floor;
    if (y > half) {
      break;
    }
    const double x2 = x * x;
    s44 += x2 * x2;
    s45 += x2 * x2 * x;
    s55 += x2 * x2 * x2;
    sy2 += y * x2;
    sy3 += y * x2 * x;
    ++used;
  }
  if (used == 0) {
    // The next image already stands above half the barrier: the parabola
    // through it is all the band says about this well.
    const size_t i = leftEnd ? 1 : n - 2;
    const double x = leftEnd ? s[i] - s.front() : s.back() - s[i];
    return 2.0 * (v[i] - floor) / (x * x);
  }
  if (used == 1) {
    return 2.0 * sy2 / s44;
  }
  const double det = s44 * s55 - s45 * s45;
  const double a = (sy2 * s55 - sy3 * s45) / det;
  return 2.0 * a;
}

wkbAction()

double eonc::tunneling::wkbAction	(	const Profile &	p,
		double	energy,
		int	points )

(1/hbar) integral sqrt(2 (V(s) - E)) ds over the path where V > E.

Definition at line 151 of file Tunneling.cpp.

                                                              {
  const double a = p.s().front();
  const double b = p.s().back();
  const double h = (b - a) / (points - 1);
  double sum = 0.0;
  for (int i = 0; i < points; ++i) {
    const double gap = p(a + i * h) - energy;
    const double f = gap > 0.0 ? std::sqrt(2.0 * gap) : 0.0;
    sum += (i == 0 || i == points - 1) ? 0.5 * f : f;
  }
  return sum * h / kHbar;
}

wkbLogRateAlongPath()

double eonc::tunneling::wkbLogRateAlongPath	(	const Profile &	profile,
		double	beta,
		double	hwReactant )

ln(k), k in 1/time, for the one-dimensional thermal rate along profile.

hwReactant is hbar omega of the reactant well, in eV. Below the barrier the transmission is the WKB factor; above it, the parabolic continuation.

Definition at line 386 of file Tunneling.cpp.

                                              {
  if (!(beta > 0.0) || !(hwReactant > 0.0)) {
    throw std::invalid_argument(
        "wkbLogRateAlongPath: beta and hbar omega must be positive");
  }
  const double vReactant = profile.v().front();
  const double s0 = profile.s().front();
  const double s1 = profile.s().back();
  double vTop = vReactant;
  double sTop = s0;
  const int grid = 2000;
  for (int k = 0; k <= grid; ++k) {
    const double sk = s0 + (s1 - s0) * static_cast<double>(k) / grid;
    const double vk = profile(sk);
    if (vk >= vTop) {
      vTop = vk;
      sTop = sk;
    }
  }
  const double barrier = vTop - vReactant;
  if (!(barrier > 0.0)) {
    throw std::invalid_argument(
        "wkbLogRateAlongPath: no barrier above the reactant");
  }
  double ds = 1e-3 * (s1 - s0);
  ds = std::min(ds, std::min(sTop - s0, s1 - sTop));
  if (!(ds > 0.0)) {
    throw std::invalid_argument(
        "wkbLogRateAlongPath: the barrier top is at an end of the path");
  }
  const double curvature =
      std::max(1e-12, -(profile(sTop + ds) - 2.0 * vTop + profile(sTop - ds)) /
                          (ds * ds));
  const double hwBarrier = kHbar * std::sqrt(curvature);
  // int P(E) exp(-beta E) dE from the reactant up to where the Boltzmann
  // factor has died. E is measured from the reactant.
  const double eMax = barrier + 40.0 / beta;
  const int points = 600;
  auto logAdd = [](double a, double b) {
    if (a == -std::numeric_limits<double>::infinity()) {
      return b;
    }
    const double m = std::max(a, b);
    return m + std::log(std::exp(a - m) + std::exp(b - m));
  };
  double logTerms = -std::numeric_limits<double>::infinity();
  double prevLog = -std::numeric_limits<double>::infinity();
  double prevE = 0.0;
  for (int k = 0; k <= points; ++k) {
    const double energy = eMax * static_cast<double>(k) / points;
    double theta = 0.0;
    if (energy < barrier) {
      theta = wkbAction(profile, vReactant + energy);
    } else {
      theta = -std::numbers::pi * (energy - barrier) / hwBarrier;
    }
    const double logP =
        theta > 20.0 ? -2.0 * theta : -std::log1p(std::exp(2.0 * theta));
    const double logF = logP - beta * energy;
    if (k > 0) {
      const double segment =
          std::log(0.5 * (energy - prevE)) + logAdd(prevLog, logF);
      logTerms = logAdd(logTerms, segment);
    }
    prevLog = logF;
    prevE = energy;
  }
  const double logFlux = logTerms - std::log(2.0 * std::numbers::pi * kHbar);
  return logFlux + std::log(2.0 * std::sinh(0.5 * beta * hwReactant));
}

wkbSplitting()

Splitting eonc::tunneling::wkbSplitting	(	const Profile &	p,
		double	hwReactant,
		double	hwProduct )

delta0 = (hbar omega / pi) exp(-S) with the Landau and Lifshitz prefactor.

The level is the higher of the two harmonic ground states, V_i + hbar omega_i / 2, and omega is the geometric mean of the wells; for a symmetric double well both reduce to the textbook formula.

Definition at line 166 of file Tunneling.cpp.

                                                                              {
  const auto &v = p.v();
  Splitting out;
  const double top = *std::max_element(v.begin(), v.end());
  out.delta = v.back() - v.front();
  out.barrier = top - v.front();
  out.hwReactant = hwReactant;
  out.hwProduct = hwProduct;
  out.referenceEnergy =
      std::max(v.front() + 0.5 * hwReactant, v.back() + 0.5 * hwProduct);
  out.action = wkbAction(p, out.referenceEnergy);
  const double hw = std::sqrt(hwReactant * hwProduct);
  out.delta0 = hw / std::numbers::pi * std::exp(-out.action);
  out.deepWells =
      (top - v.front()) > hwReactant && (top - v.back()) > hwProduct;
  return out;
}

Variable Documentation

kBoltzmann

double eonc::tunneling::kBoltzmann = 8.617333262145177e-5

inlineconstexpr

Boltzmann constant in eV / K, correctly rounded from the exact 1.380649e-23 J / K.

Definition at line 41 of file Tunneling.h.

kHbar

double eonc::tunneling::kHbar = 0.06465415129579072

inlineconstexpr

hbar in eV^0.5 amu^0.5 Angstrom, correctly rounded from the exact SI values (h = 6.62607015e-34 J s, e = 1.602176634e-19 C) and the CODATA 2022 dalton, 1.66053906892e-27 kg.

Definition at line 37 of file Tunneling.h.

kTimeUnitSeconds

double eonc::tunneling::kTimeUnitSeconds = 1.0180505717871193e-14

inlineconstexpr

One unit of time, sqrt(amu Angstrom^2 / eV), in seconds.

Definition at line 224 of file Tunneling.h.