eonc::piqtst Namespace Reference

Classes
struct	Coordinate
	The coordinate s = n . More...
struct	Plane
struct	Rate
struct	Recrossing
struct	RecrossingOptions
struct	ScanOptions

Functions
void	integrate (std::vector< Plane > &planes)
	Trapezoid integral of the mean forces, F(s_0) = 0, with errors from independent planes.
std::vector< Plane >	scan (Potential &pot, const Coordinate &coordinate, const ScanOptions &options)
	Samples one ring per plane and integrates the mean force.
Rate	rate (const std::vector< Plane > &planes, double beta)
	k = (1/2) sqrt(2 / (pi beta)) exp(-beta F(s*)) / int_{s_0}^{s*} exp(-beta F(s)) ds, the mass-weighted coordinate carrying unit mass, beta in 1 / eV.
Recrossing	recrossing (Potential &pot, const Coordinate &coordinate, const RecrossingOptions &options)
	Bennett-Chandler transmission at s*: parents sampled with the centroid held on the plane, children with Maxwell-Boltzmann ring momenta at beta / P and the plane released, propagated by RPMD.
void	validateOptions (const instanton_options_t &o)
	Throws std::invalid_argument on an inconsistent [Instanton] pi_* key.
std::vector< std::string >	runAfterInstanton (const Parameters &params, Potential &pot, const Matter &reactant, const Matter &saddle, const MatrixXd &hSaddle, const std::vector< VectorXd > &pathQ, const std::vector< double > &temperatures, std::vector< std::pair< std::string, double > > &extras)
	[Instanton] mode rate with pi_planes > 0: the planes and rate at each temperature, and with pi_recrossing_parents > 0 the transmission factor on the top plane and k_RPMD = kappa k_PI-QTST.

Function Documentation

integrate()

void eonc::piqtst::integrate

(

std::vector< Plane > &

planes

)

Trapezoid integral of the mean forces, F(s_0) = 0, with errors from independent planes.

Fills freeEnergy and freeEnergyError.

Definition at line 86 of file PIQTST.cpp.

                                         {
  if (planes.empty()) {
    return;
  }
  const MatrixXd w = trapezoidWeights(planes);
  VectorXd f(static_cast<long>(planes.size()));
  for (size_t i = 0; i < planes.size(); ++i) {
    f(static_cast<long>(i)) = planes[i].meanForce;
  }
  const VectorXd errors = forceErrors(planes);
  const VectorXd values = w * f;
  for (size_t j = 0; j < planes.size(); ++j) {
    planes[j].freeEnergy = values(static_cast<long>(j));
    planes[j].freeEnergyError =
        propagated(w.row(static_cast<long>(j)).transpose(), errors);
  }
}

rate()

Rate eonc::piqtst::rate	(	const std::vector< Plane > &	planes,
		double	beta )

k = (1/2) sqrt(2 / (pi beta)) exp(-beta F(s*)) / int_{s_0}^{s*} exp(-beta F(s)) ds, the mass-weighted coordinate carrying unit mass, beta in 1 / eV.

The integral is the trapezoid rule over the planes; the last is s*.

Definition at line 213 of file PIQTST.cpp.

                                                       {
  const long n = static_cast<long>(planes.size());
  if (n < 2) {
    throw std::invalid_argument("piqtst: the rate needs two planes");
  }
  if (!(beta > 0.0)) {
    throw std::invalid_argument("piqtst: the rate needs a positive beta");
  }
  const MatrixXd w = trapezoidWeights(planes);
  const VectorXd errors = forceErrors(planes);
  Rate r;
  double fMin = std::numeric_limits<double>::infinity();
  for (long j = 0; j < n; ++j) {
    const double f = planes[static_cast<size_t>(j)].freeEnergy;
    if (f < fMin) {
      fMin = f;
      r.reactant = j;
    }
  }
  const long top = n - 1;
  const double fTop = planes[static_cast<size_t>(top)].freeEnergy;
  r.barrier = fTop - fMin;
  r.barrierError =
      propagated((w.row(top) - w.row(r.reactant)).transpose(), errors);
  r.firstPlaneHeight = beta * (planes.front().freeEnergy - fMin);
 
  // Z = int exp(-beta F) ds by the trapezoid rule, referenced to fMin.
  VectorXd weight(n);
  double z = 0.0;
  for (long j = 0; j < n; ++j) {
    const double left = j > 0 ? planes[static_cast<size_t>(j)].s -
                                    planes[static_cast<size_t>(j - 1)].s
                              : 0.0;
    const double right = j + 1 < n ? planes[static_cast<size_t>(j + 1)].s -
                                         planes[static_cast<size_t>(j)].s
                                   : 0.0;
    weight(j) =
        0.5 * (left + right) *
        std::exp(-beta * (planes[static_cast<size_t>(j)].freeEnergy - fMin));
    z += weight(j);
  }
  weight /= z;
  r.logRate = std::log(0.5 * std::sqrt(2.0 / (std::numbers::pi * beta))) -
              beta * (fTop - fMin) - std::log(z);
  const VectorXd gradient =
      -beta * w.row(top).transpose() + beta * (w.transpose() * weight);
  r.logRateError = propagated(gradient, errors);
  return r;
}

recrossing()

Recrossing eonc::piqtst::recrossing	(	Potential &	pot,
		const Coordinate &	coordinate,
		const RecrossingOptions &	options )

Bennett-Chandler transmission at s*: parents sampled with the centroid held on the plane, children with Maxwell-Boltzmann ring momenta at beta / P and the plane released, propagated by RPMD.

The centroid velocity along s is sdot = n . M^(1/2) v_centroid.

Definition at line 263 of file PIQTST.cpp.

                                                  {
  const long dof = 3 * c.atoms;
  const Axes ax = axes(c);
  if (o.parents < 2 || o.children < 1 || o.spacing < 1 || o.steps < 4 ||
      o.equilibration < 0) {
    throw std::invalid_argument(
        "piqtst: recrossing needs two parents, one child, a positive "
        "spacing and four steps");
  }
  const VectorXd origin = c.reference + o.s * ax.b;
  VectorXd start = origin;
  if (o.seed) {
    start = o.seed(o.s);
    if (start.size() != dof) {
      throw std::invalid_argument("piqtst: a seed has the wrong length");
    }
  }
  pathintegral::RingPolymer parent(c.atoms, c.masses, c.numbers, c.free,
                                   o.ring);
  parent.setAllBeads(start.data());
  parent.setHyperplane(ax.a, origin);
  for (long step = 0; step < o.equilibration; ++step) {
    parent.step(pot, c.box, false);
  }
 
  // The children draw every normal mode from the free-ring Boltzmann
  // distribution at beta / P, which PILE's initial momenta are, and carry
  // their own stream so parent and child noise stay independent.
  pathintegral::Options childOptions = o.ring;
  childOptions.thermostat = pathintegral::Thermostat::Pile;
  childOptions.gleFile.clear();
  childOptions.seed = o.ring.seed + 0x9E3779B97F4A7C15ULL;
  pathintegral::RingPolymer child(c.atoms, c.masses, c.numbers, c.free,
                                  childOptions);
 
  const long n = o.steps + 1;
  // Per parent: sum over children of sdot(0) h(s(t) - s*) at each time,
  // and of sdot(0) h(sdot(0)).
  std::vector<VectorXd> numerator(static_cast<size_t>(o.parents),
                                  VectorXd::Zero(n));
  std::vector<double> denominator(static_cast<size_t>(o.parents), 0.0);
  Recrossing out;
  for (long p = 0; p < o.parents; ++p) {
    for (long step = 0; step < o.spacing; ++step) {
      parent.step(pot, c.box, false);
    }
    const std::vector<VectorXd> beads = parent.beads();
    VectorXd &num = numerator[static_cast<size_t>(p)];
    double &den = denominator[static_cast<size_t>(p)];
    for (long k = 0; k < o.children; ++k) {
      child.setBeads(beads);
      child.thermalMomenta();
      std::vector<VectorXd> reversed = child.momenta();
      const double forward = ax.a.dot(child.centroidVelocity());
      for (int sign = 0; sign < 2; ++sign) {
        if (sign == 1) {
          for (auto &v : reversed) {
            v = -v;
          }
          child.setBeads(beads);
          child.setMomenta(reversed);
        }
        const double sdot = sign == 0 ? forward : -forward;
        const double flux = sdot > 0.0 ? sdot : 0.0;
        den += flux;
        num(0) += flux;
        for (long i = 1; i < n; ++i) {
          child.nveStep(pot, c.box);
          const double s = ax.a.dot(child.centroid() - c.reference);
          if (s > o.s) {
            num(i) += sdot;
          }
        }
        ++out.trajectories;
      }
    }
  }
 
  VectorXd total = VectorXd::Zero(n);
  double totalDen = 0.0;
  for (long p = 0; p < o.parents; ++p) {
    total += numerator[static_cast<size_t>(p)];
    totalDen += denominator[static_cast<size_t>(p)];
  }
  if (!(totalDen > 0.0)) {
    throw std::runtime_error("piqtst: no child left the plane forward");
  }
  out.time.resize(static_cast<size_t>(n));
  out.kappa.resize(static_cast<size_t>(n));
  for (long i = 0; i < n; ++i) {
    out.time[static_cast<size_t>(i)] = static_cast<double>(i) * o.ring.dt;
    out.kappa[static_cast<size_t>(i)] = total(i) / totalDen;
  }
 
  // Plateau numerators per parent, averaged over the last quarter, and the
  // leave-one-parent-out ratios.
  const long first = n - 1 - o.steps / 4;
  const double span = static_cast<double>(n - first);
  std::vector<double> plateauNum(static_cast<size_t>(o.parents), 0.0);
  double sumNum = 0.0;
  for (long p = 0; p < o.parents; ++p) {
    plateauNum[static_cast<size_t>(p)] =
        numerator[static_cast<size_t>(p)].tail(n - first).sum() / span;
    sumNum += plateauNum[static_cast<size_t>(p)];
  }
  out.plateau = sumNum / totalDen;
  const double np = static_cast<double>(o.parents);
  std::vector<double> leaveOut(static_cast<size_t>(o.parents), 0.0);
  double meanLeave = 0.0;
  for (long p = 0; p < o.parents; ++p) {
    leaveOut[static_cast<size_t>(p)] =
        (sumNum - plateauNum[static_cast<size_t>(p)]) /
        (totalDen - denominator[static_cast<size_t>(p)]);
    meanLeave += leaveOut[static_cast<size_t>(p)];
  }
  meanLeave /= np;
  double var = 0.0;
  for (const double v : leaveOut) {
    var += (v - meanLeave) * (v - meanLeave);
  }
  out.plateauError = std::sqrt((np - 1.0) / np * var);
  out.batches = parent.batches() + child.batches();
  return out;
}

runAfterInstanton()

std::vector< std::string > eonc::piqtst::runAfterInstanton	(	const Parameters &	params,
		Potential &	pot,
		const Matter &	reactant,
		const Matter &	saddle,
		const MatrixXd &	hSaddle,
		const std::vector< VectorXd > &	pathQ,
		const std::vector< double > &	temperatures,
		std::vector< std::pair< std::string, double > > &	extras )

[Instanton] mode rate with pi_planes > 0: the planes and rate at each temperature, and with pi_recrossing_parents > 0 the transmission factor on the top plane and k_RPMD = kappa k_PI-QTST.

hSaddle and the band in pathQ are in the instanton's mass-weighted coordinates over the free atoms, measured from the reactant; saddle is aligned to the reactant. Appends results.dat keys for the last (lowest) temperature to extras and returns the files written.

Definition at line 167 of file PIQTSTJob.cpp.

                                                                 {
  const auto &o = params.instanton_options();
  validateOptions(o);
  const double timeUnit = params.constants().timeUnit;
  const Embedding emb(reactant);
  const VectorXd qSaddle = emb.toQ(reactant, saddle);
  const double qNorm = qSaddle.norm();
  if (!(qNorm > 0.0)) {
    throw std::runtime_error("piqtst: the saddle sits on the reactant");
  }
  const VectorXd lineDir = qSaddle / qNorm;
 
  // The plane normal. The saddle's unstable mode is the dividing surface
  // the rate is evaluated on; one normal for every plane keeps s a linear
  // coordinate, so the mean force integrates to its free energy.
  VectorXd dir = lineDir;
  std::string chosen = "line";
  if (o.pi_direction == "mode") {
    const Eigen::SelfAdjointEigenSolver<MatrixXd> es(
        0.5 * (hSaddle + hSaddle.transpose()));
    if (es.info() == Eigen::Success && es.eigenvalues()(0) < 0.0) {
      VectorXd v = es.eigenvectors().col(0);
      if (v.dot(lineDir) < 0.0) {
        v = -v;
      }
      if (v.dot(lineDir) >= 0.5) {
        dir = v;
        chosen = "mode";
      } else {
        EONC_LOG_WARNING("[Instanton] piqtst: the unstable mode makes {:.3f} "
                         "rad with the reactant-saddle line; planes follow "
                         "the line",
                         std::acos(v.dot(lineDir)));
      }
    } else {
      EONC_LOG_WARNING("[Instanton] piqtst: the saddle Hessian has no "
                       "negative eigenvalue; planes follow the line");
    }
  }
  if (emb.freeAtoms.size() == static_cast<size_t>(reactant.numberOfAtoms()) &&
      !reactant.getPeriodic()) {
    EONC_LOG_WARNING("[Instanton] piqtst: no atom is fixed; the planes are "
                     "fixed in the reactant's frame and a rotation of the "
                     "cluster moves the centroid between them");
  }
  const double sStar = dir.dot(qSaddle);
  const long nPlanes = o.pi_planes;
  const double s0 = -o.pi_reactant_extent * sStar;
  std::vector<double> sPlanes;
  for (long j = 0; j < nPlanes; ++j) {
    sPlanes.push_back(s0 + (sStar - s0) * static_cast<double>(j) /
                               static_cast<double>(nPlanes - 1));
  }
 
  Coordinate c;
  c.atoms = reactant.numberOfAtoms();
  c.masses.resize(static_cast<size_t>(c.atoms));
  c.numbers.resize(static_cast<size_t>(c.atoms));
  c.free.assign(static_cast<size_t>(3 * c.atoms), 0);
  const VectorXi z = reactant.getAtomicNrs();
  for (long i = 0; i < c.atoms; ++i) {
    c.masses[static_cast<size_t>(i)] = reactant.getMass(i);
    c.numbers[static_cast<size_t>(i)] = z(i);
  }
  for (const long i : emb.freeAtoms) {
    for (int a = 0; a < 3; ++a) {
      c.free[static_cast<size_t>(3 * i + a)] = 1;
    }
  }
  c.reference = flatten(reactant.getPositions());
  c.direction = emb.full(dir);
  const Matrix3d cell =
      reactant.getPeriodic() ? reactant.getCell() : Matrix3d::Zero().eval();
  c.box = cell.data();
 
  // Seeds: the band where it crosses the plane, else the straight line
  // through the reactant and the saddle.
  std::vector<double> bandS;
  for (const auto &q : pathQ) {
    bandS.push_back(dir.dot(q));
  }
  auto seed = [&](double s) -> VectorXd {
    for (size_t k = 0; k + 1 < pathQ.size(); ++k) {
      const double lo = bandS[k];
      const double hi = bandS[k + 1];
      if ((lo <= s && s <= hi) || (hi <= s && s <= lo)) {
        const double t = hi != lo ? (s - lo) / (hi - lo) : 0.0;
        return emb.cartesian(c.reference,
                             (1.0 - t) * pathQ[k] + t * pathQ[k + 1]);
      }
    }
    return emb.cartesian(c.reference, (s / sStar) * qSaddle);
  };
 
  const double vReactant = Matter(reactant).getPotentialEnergy();
  const double vSaddle = Matter(saddle).getPotentialEnergy();
  double effective = std::numeric_limits<double>::quiet_NaN();
  for (const auto &kv : extras) {
    if (kv.first == "barrier_effective_instanton") {
      effective = kv.second;
    }
  }
  EONC_LOG_INFO("[Instanton] piqtst: {} planes from s = {:.4f} to s* = {:.4f} "
                "amu^0.5 A along the {}, {} beads, {} + {} steps of {:.4g} fs",
                nPlanes, s0, sStar, chosen == "mode" ? "unstable mode" : "line",
                o.pi_beads, o.pi_equilibration_steps, o.pi_sampling_steps,
                o.pi_time_step);
 
  const std::string tableFile = "rate_piqtst.dat";
  std::ofstream table(tableFile);
  if (!table) {
    throw std::runtime_error("piqtst: cannot write " + tableFile);
  }
  const bool kappaOn = o.pi_recrossing_parents > 0;
  table << "# T_K s_amu05A dF_ds_eV_per_amu05A dF_ds_error F_eV F_error_eV "
           "spread_max_A";
  if (kappaOn) {
    table << " kappa kappa_error ln_k_rpmd_s ln_k_rpmd_s_error";
  }
  table << '\n';
  table << std::setprecision(10);
  std::vector<std::string> files{tableFile};
 
  std::vector<double> sorted = temperatures;
  std::sort(sorted.begin(), sorted.end(), std::greater<>());
  const double logSecond = std::log(tunneling::kTimeUnitSeconds);
  for (size_t ti = 0; ti < sorted.size(); ++ti) {
    const double t = sorted[ti];
    const bool last = ti + 1 == sorted.size();
    const double beta = 1.0 / (tunneling::kBoltzmann * t);
    ScanOptions so;
    so.planes = sPlanes;
    so.equilibration = o.pi_equilibration_steps;
    so.production = o.pi_sampling_steps;
    so.blocks = 10;
    so.ring.beads = o.pi_beads;
    so.ring.temperature = t;
    so.ring.kB = tunneling::kBoltzmann;
    so.ring.hbar = tunneling::kHbar;
    so.ring.dt = o.pi_time_step / timeUnit;
    so.ring.thermostat = o.pi_thermostat == "piglet"
                             ? pathintegral::Thermostat::Piglet
                             : pathintegral::Thermostat::Pile;
    so.ring.pileTau = o.pi_pile_tau / timeUnit;
    so.ring.pileScale = o.pi_pile_scale;
    so.ring.gleFile = o.pi_gle_file;
    so.ring.seed = static_cast<std::uint64_t>(o.pi_seed);
    so.seed = seed;
    const std::vector<Plane> planes = scan(pot, c, so);
    const Rate r = rate(planes, beta);
    const double lnPerSecond = r.logRate - logSecond;
 
    Recrossing kappa;
    double lnRpmd = std::numeric_limits<double>::quiet_NaN();
    double lnRpmdError = std::numeric_limits<double>::quiet_NaN();
    if (kappaOn) {
      RecrossingOptions ro;
      ro.s = planes.back().s;
      ro.equilibration = o.pi_equilibration_steps;
      ro.parents = o.pi_recrossing_parents;
      ro.spacing = o.pi_recrossing_spacing;
      ro.children = o.pi_recrossing_children;
      ro.steps = std::lround(o.pi_recrossing_time / o.pi_time_step);
      ro.ring = so.ring;
      ro.seed = seed;
      kappa = recrossing(pot, c, ro);
      if (kappa.plateau > 0.0) {
        lnRpmd = lnPerSecond + std::log(kappa.plateau);
        lnRpmdError =
            std::hypot(r.logRateError, kappa.plateauError / kappa.plateau);
      } else {
        EONC_LOG_WARNING("[Instanton] piqtst {:.4g} K: the transmission "
                         "factor {:.4g} +- {:.2g} is not positive; ln k_RPMD "
                         "is undefined (raise pi_recrossing_parents)",
                         t, kappa.plateau, kappa.plateauError);
      }
      EONC_LOG_INFO("[Instanton] piqtst {:.4g} K: transmission factor "
                    "kappa = {:.4f} +- {:.4f} from {} trajectories, "
                    "ln(k_RPMD s) = {:.4f} +- {:.4f}",
                    t, kappa.plateau, kappa.plateauError, kappa.trajectories,
                    lnRpmd, lnRpmdError);
      std::vector<std::string> curveFiles;
      if (last) {
        curveFiles.push_back("kappa_piqtst.dat");
      }
      if (sorted.size() > 1) {
        curveFiles.push_back("kappa_piqtst_" + kelvinTag(t) + ".dat");
      }
      for (const auto &file : curveFiles) {
        std::ofstream curve(file);
        if (!curve) {
          throw std::runtime_error("piqtst: cannot write " + file);
        }
        curve << "# t_fs kappa\n" << std::setprecision(10);
        for (size_t i = 0; i < kappa.time.size(); ++i) {
          curve << kappa.time[i] * timeUnit << ' ' << kappa.kappa[i] << '\n';
        }
        files.push_back(file);
      }
    }
 
    std::vector<std::string> conFiles;
    if (last) {
      conFiles.push_back("piqtst_planes.con");
    }
    if (sorted.size() > 1) {
      conFiles.push_back("piqtst_planes_" + kelvinTag(t) + ".con");
    }
    for (size_t j = 0; j < planes.size(); ++j) {
      const Plane &p = planes[j];
      double largest = 0.0;
      for (const double sp : p.spread) {
        largest = std::max(largest, sp);
      }
      table << t << ' ' << p.s << ' ' << p.meanForce << ' ' << p.meanForceError
            << ' ' << p.freeEnergy << ' ' << p.freeEnergyError << ' '
            << largest;
      if (kappaOn) {
        table << ' ' << kappa.plateau << ' ' << kappa.plateauError << ' '
              << lnRpmd << ' ' << lnRpmdError;
      }
      table << '\n';
      Matter frame(reactant);
      AtomMatrix pos = frame.getPositions();
      for (long i = 0; i < c.atoms; ++i) {
        for (int a = 0; a < 3; ++a) {
          pos(i, a) = p.centroid(3 * i + a);
        }
      }
      frame.setPositions(pos);
      io::ConFrameMetadata meta;
      meta.frame_index = static_cast<uint64_t>(j);
      meta.write_con_forces = false;
      meta.spreads = p.spread;
      meta.scalars = {{"piqtst_temperature_K", t},
                      {"piqtst_s", p.s},
                      {"piqtst_dF_ds", p.meanForce},
                      {"piqtst_dF_ds_error", p.meanForceError},
                      {"piqtst_F", p.freeEnergy},
                      {"piqtst_F_error", p.freeEnergyError},
                      {"beads", static_cast<double>(o.pi_beads)},
                      {"spread_max", largest}};
      for (const auto &file : conFiles) {
        if (!io::io_ok(frame.matter2con(file, j > 0, &meta))) {
          throw std::runtime_error("piqtst: cannot write " + file);
        }
      }
    }
    for (const auto &file : conFiles) {
      files.push_back(file);
    }
 
    const Plane &top = planes.back();
    EONC_LOG_INFO("[Instanton] piqtst {:.4g} K: free-energy barrier {:.5f} "
                  "+- {:.5f} eV (classical {:.5f} eV, instanton effective "
                  "{:.5f} eV), ln(k s) = {:.4f} +- {:.4f}",
                  t, r.barrier, r.barrierError, vSaddle - vReactant, effective,
                  lnPerSecond, r.logRateError);
    if (r.firstPlaneHeight < 5.0) {
      EONC_LOG_WARNING("[Instanton] piqtst {:.4g} K: the first plane is "
                       "{:.3g} kT above the reactant minimum of F; the "
                       "reactant integral is cut there (raise "
                       "pi_reactant_extent)",
                       t, r.firstPlaneHeight);
    }
    if (std::abs(top.meanForce) > 3.0 * top.meanForceError) {
      EONC_LOG_WARNING("[Instanton] piqtst {:.4g} K: dF/ds at the saddle "
                       "plane is {:.4g} +- {:.2g} eV / (amu^0.5 A); the "
                       "maximum of F is not at s*",
                       t, top.meanForce, top.meanForceError);
    }
    if (!last) {
      continue;
    }
    extras.emplace_back("piqtst_temperature_K", t);
    extras.emplace_back("piqtst_planes", static_cast<double>(nPlanes));
    extras.emplace_back("piqtst_beads", static_cast<double>(o.pi_beads));
    extras.emplace_back("piqtst_s_star", sStar);
    extras.emplace_back("piqtst_dF_ds_star", top.meanForce);
    extras.emplace_back("piqtst_dF_ds_star_error", top.meanForceError);
    extras.emplace_back("piqtst_first_plane_kT", r.firstPlaneHeight);
    extras.emplace_back("barrier_piqtst", r.barrier);
    extras.emplace_back("barrier_piqtst_error", r.barrierError);
    extras.emplace_back("rate_piqtst", std::exp(lnPerSecond));
    extras.emplace_back("rate_piqtst_log", lnPerSecond);
    extras.emplace_back("rate_piqtst_log_error", r.logRateError);
    if (kappaOn) {
      extras.emplace_back("piqtst_kappa", kappa.plateau);
      extras.emplace_back("piqtst_kappa_error", kappa.plateauError);
      extras.emplace_back("rate_rpmd", std::exp(lnRpmd));
      extras.emplace_back("rate_rpmd_log", lnRpmd);
      extras.emplace_back("rate_rpmd_log_error", lnRpmdError);
    }
  }
  return files;
}

scan()

std::vector< Plane > eonc::piqtst::scan	(	Potential &	pot,
		const Coordinate &	coordinate,
		const ScanOptions &	options )

Samples one ring per plane and integrates the mean force.

One ring is carried from plane to plane: its centroid moves to the next seed and the internal modes keep their state.

Definition at line 104 of file PIQTST.cpp.

                                              {
  const long dof = 3 * c.atoms;
  const Axes ax = axes(c);
  const VectorXd &a = ax.a;
  const VectorXd &b = ax.b;
  if (o.planes.size() < 2) {
    throw std::invalid_argument("piqtst: at least two planes are needed");
  }
  for (size_t j = 1; j < o.planes.size(); ++j) {
    if (!(o.planes[j] > o.planes[j - 1])) {
      throw std::invalid_argument("piqtst: plane positions must ascend");
    }
  }
  if (o.equilibration < 0 || o.production < 2 || o.blocks < 2 ||
      o.production < o.blocks) {
    throw std::invalid_argument(
        "piqtst: sampling needs at least as many steps as blocks, and two "
        "blocks");
  }
  const double aNorm = a.norm();
  auto seedAt = [&](double s) -> VectorXd {
    if (o.seed) {
      VectorXd x = o.seed(s);
      if (x.size() != dof) {
        throw std::invalid_argument("piqtst: a seed has the wrong length");
      }
      return x;
    }
    return c.reference + s * b;
  };
 
  pathintegral::RingPolymer ring(c.atoms, c.masses, c.numbers, c.free, o.ring);
  const long beads = o.ring.beads;
  const long blockSize = o.production / o.blocks;
  std::vector<Plane> out;
  out.reserve(o.planes.size());
  for (size_t j = 0; j < o.planes.size(); ++j) {
    const double s = o.planes[j];
    const VectorXd origin = c.reference + s * b;
    const VectorXd target = seedAt(s);
    if (j == 0) {
      ring.setAllBeads(target.data());
    } else {
      const VectorXd shift = target - ring.centroid();
      std::vector<VectorXd> moved = ring.beads();
      for (auto &q : moved) {
        q += shift;
      }
      ring.setBeads(moved);
    }
    ring.setHyperplane(a, origin);
 
    const long batches0 = ring.batches();
    for (long step = 0; step < o.equilibration; ++step) {
      ring.step(pot, c.box, false);
    }
    Plane plane;
    plane.s = s;
    plane.centroid = VectorXd::Zero(dof);
    VectorXd spread2 = VectorXd::Zero(dof);
    double sum = 0.0;
    double blockSum = 0.0;
    std::vector<double> blockMeans;
    for (long step = 0; step < o.production; ++step) {
      ring.resetAverages();
      ring.step(pot, c.box, true);
      const double fn = ring.meanForce();
      sum += fn;
      blockSum += fn;
      if ((step + 1) % blockSize == 0 &&
          static_cast<long>(blockMeans.size()) < o.blocks) {
        blockMeans.push_back(blockSum / static_cast<double>(blockSize));
        blockSum = 0.0;
      }
      const VectorXd centroid = ring.centroid();
      plane.centroid += centroid;
      for (const auto &q : ring.beads()) {
        spread2 += (q - centroid).cwiseAbs2();
      }
    }
    const double steps = static_cast<double>(o.production);
    plane.centroid /= steps;
    spread2 /= steps * static_cast<double>(beads);
    plane.spread.resize(static_cast<size_t>(dof));
    for (long i = 0; i < dof; ++i) {
      plane.spread[static_cast<size_t>(i)] = std::sqrt(spread2(i));
    }
    const double mean = sum / steps;
    double blockMean = 0.0;
    for (const double m : blockMeans) {
      blockMean += m;
    }
    const double nb = static_cast<double>(blockMeans.size());
    blockMean /= nb;
    double var = 0.0;
    for (const double m : blockMeans) {
      var += (m - blockMean) * (m - blockMean);
    }
    var /= nb * (nb - 1.0);
    plane.meanForce = -mean / aNorm;
    plane.meanForceError = std::sqrt(var) / aNorm;
    plane.batches = ring.batches() - batches0;
    out.push_back(std::move(plane));
  }
  integrate(out);
  return out;
}

validateOptions()

void eonc::piqtst::validateOptions

(

const instanton_options_t &

o

)

Throws std::invalid_argument on an inconsistent [Instanton] pi_* key.

Definition at line 31 of file PIQTSTJob.cpp.

                                                   {
  if (o.pi_planes < 0) {
    throw std::invalid_argument("[Instanton] pi_planes must not be negative");
  }
  if (o.pi_planes == 0) {
    return;
  }
  if (o.mode != "rate") {
    throw std::invalid_argument("[Instanton] pi_planes needs mode = rate");
  }
  if (o.pi_planes < 2) {
    throw std::invalid_argument(
        "[Instanton] pi_planes must be 0 or at least 2");
  }
  if (o.pi_beads < 1) {
    throw std::invalid_argument("[Instanton] pi_beads must be positive");
  }
  if (o.pi_equilibration_steps < 0 || o.pi_sampling_steps < 20) {
    throw std::invalid_argument(
        "[Instanton] pi_equilibration_steps must not be negative and "
        "pi_sampling_steps must be at least 20 (ten blocks of two)");
  }
  if (!(o.pi_time_step > 0.0) || !(o.pi_pile_tau > 0.0) ||
      !(o.pi_pile_scale > 0.0)) {
    throw std::invalid_argument("[Instanton] pi_time_step, pi_pile_tau and "
                                "pi_pile_scale must be positive");
  }
  if (o.pi_seed < 0) {
    throw std::invalid_argument("[Instanton] pi_seed must not be negative");
  }
  if (o.pi_thermostat != "pile" && o.pi_thermostat != "piglet") {
    throw std::invalid_argument(
        "[Instanton] pi_thermostat must be pile or piglet, not " +
        o.pi_thermostat);
  }
  if (o.pi_thermostat == "piglet" && o.pi_beads > 1 && o.pi_gle_file.empty()) {
    throw std::invalid_argument(
        "[Instanton] pi_thermostat = piglet needs pi_gle_file");
  }
  if (o.pi_direction != "mode" && o.pi_direction != "line") {
    throw std::invalid_argument(
        "[Instanton] pi_direction must be mode or line, not " + o.pi_direction);
  }
  if (!(o.pi_reactant_extent >= 0.0)) {
    throw std::invalid_argument(
        "[Instanton] pi_reactant_extent must not be negative");
  }
  if (o.pi_recrossing_parents < 0 || o.pi_recrossing_parents == 1) {
    throw std::invalid_argument(
        "[Instanton] pi_recrossing_parents must be 0 or at least 2");
  }
  if (o.pi_recrossing_parents == 0) {
    return;
  }
  if (o.pi_recrossing_children < 1 || o.pi_recrossing_spacing < 1) {
    throw std::invalid_argument("[Instanton] pi_recrossing_children and "
                                "pi_recrossing_spacing must be positive");
  }
  if (!(o.pi_recrossing_time >= 4.0 * o.pi_time_step)) {
    throw std::invalid_argument(
        "[Instanton] pi_recrossing_time must be at least four pi_time_step");
  }
}