Tunneling.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
 
// One-dimensional tunnelling along a mass-weighted minimum energy path.
//
// A two-level system in a glass is a pair of adjacent minima that the system
// tunnels between at about one kelvin. Its splitting follows from a WKB
// integral along the mass-weighted path joining the minima (Damart and
// Rodney; Khomenko et al.; Mocanu et al.). Energies are in eV, lengths in
// Angstrom and masses in amu, so mass-weighted lengths are in
// amu^0.5 Angstrom and hbar is kHbar.
 
#include "Matter.h"
 
#include <array>
#include <functional>
#include <limits>
#include <memory>
#include <string>
#include <vector>
 
namespace eonc::tunneling {

inline constexpr double kHbar = 0.06465415129579072;

inline constexpr double kBoltzmann = 8.617333262145177e-5;

double massWeightedDistance(const Matter &a, const Matter &b);

std::vector<double>
massWeightedPath(const std::vector<std::shared_ptr<Matter>> &band);

class Profile {
public:
  Profile(std::vector<double> s, std::vector<double> v);
  double operator()(double x) const;
  const std::vector<double> &s() const { return s_; }
  const std::vector<double> &v() const { return v_; }
 
private:
  std::vector<double> s_, v_, m_;
};

double wellCurvature(const Profile &p, bool leftEnd);

double hbarOmega(double curvature);

double wkbAction(const Profile &p, double energy, int points = 4001);
 
struct Splitting {
  double delta = 0.0;           
  double barrier = 0.0;         
  double hwReactant = 0.0;      
  double hwProduct = 0.0;       
  double referenceEnergy = 0.0; 
  double action = 0.0;          
  double delta0 = 0.0;          
  bool deepWells = false;
  double tlsEnergy() const; 
};

Splitting wkbSplitting(const Profile &p, double hwReactant, double hwProduct);

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

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

double wkbLogRateAlongPath(const Profile &profile, double beta,
                           double hwReactant);
 
// Ring-polymer instanton for the splitting between two minima.
//
// A path of P + 1 beads in mass-weighted coordinates q runs from one minimum
// to the other in imaginary time beta hbar, the ends fixed at the minima. The
// instanton is the path that minimises the discretised Euclidean action
//
//   S = sum_j |q_{j+1} - q_j|^2 / (2 dtau) + dtau sum_j' V(q_j),
//
// dtau = beta hbar / P, with half weight on the two ends (trapezoid). The
// splitting follows from the ratio of the off-diagonal to the diagonal
// imaginary-time propagator, both in the same steepest-descent
// approximation:
//
//   delta0 = 2 hbar sqrt(S0 / (2 pi hbar dtau))
//            sqrt(det J_well / det' J) exp(-(S - S_well) / hbar),
//
// with J the Hessian of S over the interior beads, det' leaving out its
// zero mode (the kink's position in imaginary time), S0 the integral of
// |dq/dtau|^2 and J_well the same Hessian with every bead at a minimum.
// Time runs in amu^0.5 Angstrom eV^-0.5, so hbar is kHbar. Where the minimum
// energy path curves, the instanton cuts the corner and follows the
// transverse zero-point energy, which a one-dimensional WKB integral along
// the path cannot.
 
struct InstantonOptions {
  long beads = 256;             
  double betaHbarOmega = 30.0;  
  long maxIterations = 5000;    
  double forceTolerance = 1e-4; 
  long memory = 20;             
};

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

using BeadHessian = std::function<MatrixXd(long j, const VectorXd &q)>;
 
struct Instanton {
  std::vector<VectorXd> path;   
  std::vector<double> energies; 
  double betaHbar = 0.0;        
  double dtau = 0.0;            
  double action = 0.0;          
  double s0 = 0.0;              
  double zeroMode = 0.0;        
  double delta0 = 0.0;          
  double asymmetry = 0.0;       
  long iterations = 0;
  bool converged = false;
  bool symmetricEnough = false;
  double modeSeparation = 0.0;
};

double pathOmega(const MatrixXd &hessStart, const MatrixXd &hessEnd,
                 const VectorXd &start, const VectorXd &end);

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

void instantonSplitting(Instanton &inst, const BeadHessian &hessian,
                        const MatrixXd &hessStart, const MatrixXd &hessEnd);
 
// Ring-polymer instanton for the thermal rate below the crossover
// temperature (Richardson and Althorpe, J. Chem. Phys. 131, 214106 (2009)).
//
// A closed ring of N beads in mass-weighted coordinates q, beta_N = beta / N,
// has the potential
//
//   U_N = sum_j V(q_j) + sum_j |q_{j+1} - q_j|^2 / (2 beta_N^2 hbar^2),
//
// q_N = q_0. Below T_c = hbar omega_b / (2 pi kB), omega_b the imaginary
// frequency at the saddle, the instanton is a first-order saddle of U_N:
// one negative mode, and one zero mode that cycles the beads. The rate is
//
//   k Z_r = (1 / (beta_N hbar)) sqrt(B_N / (2 pi beta_N hbar^2))
//           prod'_k 1 / (beta_N hbar |omega_k|) exp(-beta_N U_N),
//
// B_N = sum_j |q_{j+1} - q_j|^2, omega_k^2 the eigenvalues of the
// mass-weighted Hessian of U_N with the zero mode left out, and Z_r the
// ring-polymer partition function of the harmonic reactant,
//
//   Z_r = exp(-beta V_r) prod_{k=0}^{N-1} prod_m
//         1 / (beta_N hbar sqrt(lambda_m + 4 sin^2(pi k / N) / (beta_N
//         hbar)^2)),
//
// lambda_m the eigenvalues of the reactant's mass-weighted Hessian. Above
// T_c the ring collapses onto the saddle and the formula no longer applies.

inline constexpr double kTimeUnitSeconds = 1.0180505717871193e-14;

double crossoverTemperature(const MatrixXd &hessSaddle);

double parabolicFactor(double temperature, double crossover);

double harmonicTstLogRate(const MatrixXd &hessReactant,
                          const MatrixXd &hessSaddle, double beta,
                          double barrier, long rigidModes);

double quantumHarmonicTstLogRate(const MatrixXd &hessReactant,
                                 const MatrixXd &hessSaddle, double beta,
                                 double barrier, long rigidModes);
 
struct RateInstantonOptions {
  long beads = 32;              
  long maxIterations = 1000;    
  double forceTolerance = 1e-3; 
  long lanczosFirst = 30;       
  long lanczosRestart = 6;      
  double lanczosStep = 1e-4;    
  double maxStep = 0.05;        
  long memory = 10;             
  bool halfRing = true;
  bool checkOddSector = true;
  double energyShift = 0.0; 
  long newtonLimit = std::numeric_limits<long>::max();
  std::string initialHessians = "saddle";
  std::vector<double> rigidSqrtMasses;
  VectorXd rigidReference;
  std::array<bool, 3> rigidRotations{{false, false, false}};
};

struct RingSpectrum {
  double logDetPrime = 0.0;
  long negativeModes = 0;
  double zeroEigenvalue = 0.0;
};

RingSpectrum ringSpectrum(const std::vector<MatrixXd> &beadHessians, double c,
                          const std::vector<VectorXd> &tau);
 
struct RateInstanton {
  std::vector<VectorXd> beads;     
  std::vector<double> energies;    
  double beta = 0.0;               
  double betaN = 0.0;              
  double temperature = 0.0;        
  double crossover = 0.0;          
  double ringPotential = 0.0;      
  double bN = 0.0;                 
  double negativeEigenvalue = 0.0; 
  double zeroEigenvalue = 0.0;     
  long negativeModes = 0;          
  long iterations = 0;
  bool converged = false;
  double logRateTimesZr = 0.0; 
  double logZr = 0.0;          
  double logRate = 0.0;        
  double rate = 0.0;           
  double effectiveBarrier = 0.0;
  double classicalRate = 0.0;
  double classicalLogRate = 0.0;
};

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

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

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

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

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