PIQTST.h

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** SPDX-License-Identifier: BSD-3-Clause
**
** Copyright (c) 2010--present, eOn Development Team
** All rights reserved.
**
** Repo:
** https://github.com/TheochemUI/eOn
*/
#pragma once
 
// Path-integral quantum transition-state theory (Voth, Chandler and
// Miller, J. Chem. Phys. 91, 7749 (1989)): the centroid potential of mean
// force along a linear mass-weighted coordinate s from ring polymers whose
// centroid is held on parallel planes, and the rate from its value on the
// dividing plane, optionally times the ring-polymer MD transmission factor
// on that plane (Craig and Manolopoulos, J. Chem. Phys. 122, 084106 (2005);
// Suleimanov, Allen and Green, Comput. Phys. Commun. 184, 833 (2013)).
 
#include "eon/Eigen.h"
#include "eon/ParametersOptions.h"
#include "eon/PathIntegral.h"
#include "eon/Potential.h"
 
#include <functional>
#include <string>
#include <utility>
#include <vector>
 
namespace eonc {
class Matter;
class Parameters;
} // namespace eonc
 
namespace eonc::piqtst {

struct Coordinate {
  long atoms{0};
  std::vector<double> masses; // per atom, amu
  std::vector<int> numbers;   // per atom
  std::vector<char> free;     // per coordinate, 3 * atoms
  VectorXd reference;         // Cartesian, 3 * atoms, Angstrom
  VectorXd direction;         // mass-weighted unit vector, 3 * atoms
  const double *box{nullptr}; // 3 x 3 cell, or null
};
 
struct ScanOptions {
  std::vector<double> planes;
  long equilibration{500};
  long production{2000};
  long blocks{10};
  pathintegral::Options ring;
  std::function<VectorXd(double)> seed;
};
 
struct Plane {
  double s{0.0};
  double meanForce{0.0};
  double meanForceError{0.0};
  double freeEnergy{0.0};
  double freeEnergyError{0.0};
  VectorXd centroid;
  std::vector<double> spread;
  long batches{0};
};

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

void integrate(std::vector<Plane> &planes);
 
struct Rate {
  long reactant{0};
  double barrier{0.0};
  double barrierError{0.0};
  double logRate{0.0};
  double logRateError{0.0};
  double firstPlaneHeight{0.0};
};

Rate rate(const std::vector<Plane> &planes, double beta);
 
struct RecrossingOptions {
  double s{0.0};
  long equilibration{500};
  long parents{0};
  long spacing{50};
  long children{20};
  long steps{0};
  pathintegral::Options ring;
  std::function<VectorXd(double)> seed;
};
 
struct Recrossing {
  std::vector<double> time;
  std::vector<double> kappa;
  double plateau{0.0};
  double plateauError{0.0};
  long trajectories{0};
  long batches{0};
};

Recrossing recrossing(Potential &pot, const Coordinate &coordinate,
                      const RecrossingOptions &options);

void validateOptions(const instanton_options_t &o);

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);
 
} // namespace eonc::piqtst