eonc::InstantonJob Class Reference

Ring-polymer instanton ([Instanton]): the tunnelling splitting between two minima beyond one-dimensional WKB (mode splitting), or the thermal rate through a saddle (mode rate): the ring below the crossover temperature, and the parabolic barrier factor times harmonic transition-state theory above it, written to instanton.con and results.dat. More...

#include <InstantonJob.h>

Inheritance diagram for eonc::InstantonJob:

Public Member Functions
	InstantonJob (std::unique_ptr< Parameters > parameters, Runtime &rt)
	~InstantonJob (void)=default
std::vector< std::string >	run (void) override
	Virtual run; used solely for dynamic dispatch.
Public Member Functions inherited from eonc::Job
void	adoptRuntime (std::unique_ptr< Runtime > rt)
	Take ownership of a Runtime previously passed as Runtime&.
	Job (std::unique_ptr< Parameters > parameters, Runtime &rt)
	Borrow: caller keeps Runtime alive (CLI stack / Python Session).
	Job (std::unique_ptr< Parameters > parameters, std::unique_ptr< Runtime > rt)
	Own a Runtime (one-shot makeJob / rvalue).
	Job (std::unique_ptr< Parameters > parameters)
	Own a default-constructed Runtime.
	Job (std::shared_ptr< Potential > potPassed, const Parameters &parameters)
virtual	~Job ()=default
JobType	getType ()
PotRegistry &	pots () noexcept
void	releasePotential ()
	Drop the Potential so on_destroyed is recorded before Runtime dies.

Additional Inherited Members
Protected Attributes inherited from eonc::Job
JobType	jtype
Parameters	params
std::unique_ptr< Runtime >	owned_runtime_
	Non-null when this Job owns the composition root (one-shot makeJob).
Runtime *	runtime_
	Always valid: either owned_runtime_.get() or a caller-owned Runtime.
std::shared_ptr< Potential >	pot

Detailed Description

Ring-polymer instanton ([Instanton]): the tunnelling splitting between two minima beyond one-dimensional WKB (mode splitting), or the thermal rate through a saddle (mode rate): the ring below the crossover temperature, and the parabolic barrier factor times harmonic transition-state theory above it, written to instanton.con and results.dat.

Definition at line 24 of file InstantonJob.h.

Constructor & Destructor Documentation

InstantonJob()

eonc::InstantonJob::InstantonJob ( std::unique_ptr< Parameters > parameters, Runtime & rt )

inline

Definition at line 26 of file InstantonJob.h.

      : Job(std::move(parameters), rt) {}

~InstantonJob()

eonc::InstantonJob::~InstantonJob ( void )

default

Member Function Documentation

run()

std::vector< std::string > eonc::InstantonJob::run ( void )

overridevirtual

Virtual run; used solely for dynamic dispatch.

Implements eonc::Job.

Definition at line 921 of file InstantonJob.cpp.

                                           {
  const auto &o = params.instanton_options();
  pathintegral::requireTrotterSprings(o.springs, "instanton");
  std::vector<std::string> returnFiles;
  const std::string resultsFile = "results.dat";
  const std::string pathFile = "instanton.con";
  returnFiles.push_back(resultsFile);
 
  auto reactant = std::make_unique<Matter>(pot, params);
  if (!io::io_ok(reactant->con2matter(o.reactant_filename))) {
    throw std::runtime_error("instanton: cannot read " + o.reactant_filename);
  }
  const MassWeighted mw(*reactant);
  const long n = mw.dimension();
  const VectorXd qStart = VectorXd::Zero(n);
 
  // Bead evaluations: one batch per call, through forceBatch when the
  // potential spreads a batch over calculators.
  std::vector<std::unique_ptr<Matter>> pool;
  auto evaluate = [&](const std::vector<VectorXd> &q, std::vector<double> &v,
                      std::vector<VectorXd> &grad) {
    while (pool.size() < q.size()) {
      pool.push_back(std::make_unique<Matter>(*reactant));
    }
    for (size_t j = 0; j < q.size(); ++j) {
      mw.place(q[j], *pool[j]);
    }
    v.resize(q.size());
    grad.resize(q.size());
    if (pot->supportsBatchEvaluation() && q.size() > 1) {
      const long atoms = reactant->numberOfAtoms();
      std::vector<VectorXi> nrs(q.size());
      std::vector<Matrix3d> boxes(q.size());
      std::vector<const double *> posPtr, boxPtr;
      std::vector<const int *> nrsPtr;
      std::vector<double *> frcPtr;
      for (size_t j = 0; j < q.size(); ++j) {
        nrs[j] = pool[j]->getAtomicNrs();
        boxes[j] = pool[j]->getPeriodic() ? pool[j]->getCell()
                                          : Matrix3d::Zero().eval();
      }
      for (size_t j = 0; j < q.size(); ++j) {
        posPtr.push_back(pool[j]->getPositions().data());
        nrsPtr.push_back(nrs[j].data());
        frcPtr.push_back(pool[j]->forcesData());
        boxPtr.push_back(boxes[j].data());
      }
      std::vector<double> energies(q.size()), variances(q.size());
      pot->forceBatch(static_cast<long>(q.size()), atoms, posPtr.data(),
                      nrsPtr.data(), frcPtr.data(), energies.data(),
                      variances.data(), boxPtr.data());
      for (size_t j = 0; j < q.size(); ++j) {
        pool[j]->setComputedPotential(energies[j], variances[j]);
      }
    }
    for (size_t j = 0; j < q.size(); ++j) {
      v[j] = pool[j]->getPotentialEnergy();
      grad[j] = mw.gradient(pool[j]->getForces());
    }
  };
 
  // The first call is the reactant. Its Hessian decides which rotations
  // are zero modes; later beads reuse that decision.
  std::array<bool, 3> rotationZero{{false, false, false}};
  std::array<double, 3> rotationResidual{{0.0, 0.0, 0.0}};
  bool rotationsKnown = false;
  auto hessianAt = [&](const VectorXd &q) {
    Matter m(*reactant);
    mw.place(q, m);
    Hessian h(params, &m);
    h.writeHessianFile(false);
    MatrixXd out = h.getHessian(&m, mw.freeAtoms());
    if (out.rows() != n) {
      throw std::runtime_error("instanton: a bead Hessian failed");
    }
    if (!rotationsKnown) {
      mw.markRotationZeroModes(out, mw.rigidGenerators(m), rotationZero,
                               rotationResidual);
      rotationsKnown = true;
    }
    // A finite-difference Hessian of a free structure has small nonzero
    // rigid eigenvalues of either sign; project them to zero.
    const MatrixXd rigid = mw.rigidBasis(m, rotationZero);
    if (rigid.cols() > 0) {
      const MatrixXd p = MatrixXd::Identity(n, n) - rigid * rigid.transpose();
      out = p * out * p;
    }
    return out;
  };
 
  if (o.mode == "rate") {
    return runRate(params, pot, *reactant, mw, evaluate, hessianAt,
                   rotationZero, rotationResidual);
  }
 
  auto product = std::make_unique<Matter>(pot, params);
  if (!io::io_ok(product->con2matter(o.product_filename))) {
    throw std::runtime_error("instanton: cannot read " + o.product_filename);
  }
  if (reactant->numberOfAtoms() != product->numberOfAtoms()) {
    throw std::runtime_error("instanton: the minima differ in atom count");
  }
  alignRigid(*reactant, *product);
  const VectorXd qEnd = mw.toQ(*product);
 
  // Starting path: a band from file, else the straight line.
  std::vector<VectorXd> guess;
  if (!o.initial_path.empty()) {
    const auto frames = readcon::read_all_frames(o.initial_path);
    for (const auto &frame : frames) {
      Matter m(*reactant);
      if (!io::io_ok(io::con2matter(m, frame))) {
        throw std::runtime_error("instanton: cannot read " + o.initial_path);
      }
      alignRigid(*reactant, m);
      guess.push_back(mw.toQ(m));
    }
    if (guess.size() >= 2) {
      guess.front() = qStart;
      guess.back() = qEnd;
    }
  }
 
  const MatrixXd hStart = hessianAt(qStart);
  const MatrixXd hEnd = hessianAt(qEnd);
  const double omega = tunneling::pathOmega(hStart, hEnd, qStart, qEnd);
  const double betaHbar = o.beta_hbar_omega / omega;
  tunneling::InstantonOptions opt;
  opt.beads = o.beads;
  opt.betaHbarOmega = o.beta_hbar_omega;
  opt.maxIterations = o.max_iterations;
  opt.forceTolerance = o.force_tolerance;
  EONC_LOG_INFO("[Instanton] {} beads over beta hbar = {:.4f} fs, {} degrees "
                "of freedom",
                o.beads, betaHbar * kTimeUnitFs, n);
 
  tunneling::Instanton inst = tunneling::optimizeInstanton(
      qStart, qEnd, betaHbar, guess, evaluate, opt);
  EONC_LOG_INFO("[Instanton] action {:.6f} after {} iterations{}", inst.action,
                inst.iterations, inst.converged ? "" : " (not converged)");
 
  bool splitOk = false;
  std::string failure;
  // beta |delta|: the propagator ratio reads delta0 only when the wells
  // lie within a small fraction of kB T of each other.
  const double betaAsymmetry =
      std::abs(inst.asymmetry) * betaHbar / tunneling::kHbar;
  if (inst.converged && !inst.symmetricEnough) {
    EONC_LOG_WARNING("[Instanton] beta |delta| = {:.3g}: the wells differ by "
                     "{:.4g} eV, too far for the splitting; the path and "
                     "action are written, the splitting is not",
                     betaAsymmetry, inst.asymmetry);
  }
  if (inst.converged && inst.symmetricEnough) {
    const long stride = std::max<long>(1, o.hessian_stride);
    const long P = o.beads;
    std::map<long, MatrixXd> anchors;
    auto anchor = [&](long j) -> const MatrixXd & {
      auto it = anchors.find(j);
      if (it == anchors.end()) {
        it = anchors.emplace(j, hessianAt(inst.path[static_cast<size_t>(j)]))
                 .first;
      }
      return it->second;
    };
    auto beadHessian = [&](long j, const VectorXd &) -> MatrixXd {
      if (stride == 1) {
        return anchor(j);
      }
      const long lo = 1 + ((j - 1) / stride) * stride;
      const long hi = std::min(lo + stride, P - 1);
      if (j == lo || hi == lo) {
        return anchor(lo);
      }
      const double t =
          static_cast<double>(j - lo) / static_cast<double>(hi - lo);
      return (1.0 - t) * anchor(lo) + t * anchor(hi);
    };
    try {
      tunneling::instantonSplitting(inst, beadHessian, hStart, hEnd);
      splitOk = std::isfinite(inst.delta0);
    } catch (const std::runtime_error &ex) {
      failure = ex.what();
      EONC_LOG_ERROR("[Instanton] {}", failure);
    }
  }
 
  // The path, one frame per bead; the splitting on the first frame.
  const double kelvin = tunneling::kHbar / (tunneling::kBoltzmann * betaHbar);
  Matter frame(*reactant);
  for (size_t j = 0; j < inst.path.size(); ++j) {
    mw.place(inst.path[j], frame);
    io::ConFrameMetadata meta;
    meta.frame_index = static_cast<uint64_t>(j);
    meta.energy = inst.energies[j];
    meta.write_con_forces = false;
    meta.scalars = {{"imaginary_time_fs",
                     static_cast<double>(j) * inst.dtau * kTimeUnitFs}};
    if (j == 0) {
      meta.scalars.push_back({"instanton_action", inst.action});
      meta.scalars.push_back(
          {"instanton_beta_hbar_fs", betaHbar * kTimeUnitFs});
      meta.scalars.push_back({"instanton_temperature_K", kelvin});
      meta.scalars.push_back(
          {"instanton_converged", inst.converged ? 1.0 : 0.0});
      meta.scalars.push_back({"tunnel_asymmetry", inst.asymmetry});
      if (splitOk) {
        meta.scalars.push_back({"tunnel_splitting_instanton", inst.delta0});
        meta.scalars.push_back(
            {"tls_energy_instanton", std::hypot(inst.asymmetry, inst.delta0)});
        meta.scalars.push_back({"instanton_zero_mode", inst.zeroMode});
        meta.scalars.push_back(
            {"instanton_mode_separation", inst.modeSeparation});
        meta.scalars.push_back(
            {"instanton_symmetric", inst.symmetricEnough ? 1.0 : 0.0});
      }
    }
    if (!io::io_ok(frame.matter2con(pathFile, j > 0, &meta))) {
      throw std::runtime_error("instanton: cannot write " + pathFile);
    }
  }
  returnFiles.push_back(pathFile);
  writeCentroid("instanton_centroid.con", inst.path, mw, *reactant,
                {{"instanton_temperature_K", kelvin},
                 {"instanton_converged", inst.converged ? 1.0 : 0.0}});
  returnFiles.push_back("instanton_centroid.con");
 
  // A converged path between wells too far apart is a result, not a
  // failure: the flags say why no splitting was written.
  const bool good = splitOk || (inst.converged && !inst.symmetricEnough);
  const auto status = good ? RunStatus::GOOD
                           : (inst.converged ? RunStatus::FAIL_POTENTIAL_FAILED
                                             : RunStatus::FAIL_MAX_ITERATIONS);
  auto env = JobResultEnvelope::fromMinimization(
      status, params.potential_options().potential,
      PotRegistry::get().total_force_calls(), false, 0.0);
  env.job_type = "instanton";
  env.extras.emplace_back("force_calls", static_cast<double>(env.force_calls));
  env.extras.emplace_back("instanton_iterations",
                          static_cast<double>(inst.iterations));
  env.extras.emplace_back("instanton_action", inst.action);
  env.extras.emplace_back("instanton_temperature_K", kelvin);
  env.extras.emplace_back("tunnel_asymmetry", inst.asymmetry);
  env.extras.emplace_back("instanton_beta_asymmetry", betaAsymmetry);
  env.extras.emplace_back("instanton_symmetric",
                          inst.symmetricEnough ? 1.0 : 0.0);
  if (splitOk) {
    env.extras.emplace_back("tunnel_splitting_instanton", inst.delta0);
    env.extras.emplace_back("instanton_mode_separation", inst.modeSeparation);
  }
  env.writeResultsDat(resultsFile);
  return returnFiles;
}

---

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

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