Tunnelling splittings

A two-level system (TLS) in a glass is a pair of adjacent minima that the structure tunnels between at about one kelvin. Its tunnelling splitting \(\Delta_0\) and asymmetry \(\Delta\) set the TLS energy \(E = \sqrt{\Delta^2 + \Delta_0^2}\). eOn estimates \(\Delta_0\) in two ways:

	NEB band (WKB)	job = instanton
Path	the minimum energy path	the path of least imaginary-time action
Dimensions	one, along the mass-weighted band	every free degree of freedom
Cost	the NEB itself	one batch of forces per iteration over the beads, plus a Hessian per bead
Output	first frame of neb.con	instanton.con and results.dat

With mode = rate the same job estimates the thermal rate through a saddle instead of a splitting. That path is below.

Energies are in eV, lengths in Å and masses in amu throughout, so mass-weighted lengths are in amu^0.5 Å.

WKB along a NEB band

Every nudged elastic band (NEB) job writes reaction_coordinate_mw, the mass-weighted arc length, on each frame of neb.con. The first frame also carries the band’s Wentzel-Kramers-Brillouin (WKB) estimate. The keys are:

Key	Meaning
hbar_omega_reactant, hbar_omega_product	Wells: \(\hbar\omega\) of each well along the band, from a fit of \(a s^2 + b s^3\) to the images within half the barrier
tunnel_action	Action: \(S = \hbar^{-1} \int \sqrt{2 (V(s) - E)}\, ds\) over the forbidden region
tunnel_splitting	Estimate: \(\Delta_0 = (\hbar\omega / \pi) e^{-S}\), with \(\omega\) the geometric mean of the wells
tls_energy	Energy: \(\sqrt{\Delta^2 + \Delta_0^2}\)
tunnel_deep_wells	Flag: 1 when both barriers exceed \(\hbar\omega\); below that, WKB is the wrong tool

The profile between images is a monotone cubic, so it cannot dip below the data. The level \(E\) is the higher of the two harmonic ground states. A structure without masses, or a band whose end is flat, leaves these keys out.

WKB along the band is exact in one dimension up to its semiclassical error. When the path curves, the tunnelling cuts the corner, and the transverse zero-point energy changes along the way. In the two-dimensional test valley below, both effects together put the band estimate a factor of 2.8 below the exact splitting.

The instanton job

The instanton is a path of \(P + 1\) beads in imaginary time \(\beta\hbar\), with its ends fixed at the two minima. It minimises the discretised Euclidean action

\[S = \sum_j \frac{|q_{j+1} - q_j|^2}{2\,\delta\tau} + \delta\tau \sum_j V(q_j), \qquad \delta\tau = \beta\hbar / P,\]

in mass-weighted coordinates \(q\). The splitting comes from the ratio of the off-diagonal to the diagonal imaginary-time propagator, both taken in the same steepest-descent approximation.

\[\Delta_0 = 2\hbar \sqrt{\frac{S_0}{2\pi\hbar\,\delta\tau}} \sqrt{\frac{\det J_\text{well}}{\det' J}}\; e^{-(S - S_\text{well})/\hbar}.\]

Here \(J\) is the Hessian of \(S\) over the interior beads. The prime leaves out its zero mode, the kink’s position in imaginary time. \(S_0 = \int |\dot q|^2 d\tau\), and \(J_\text{well}\) is the same Hessian with every bead at a minimum.

[Main]
job = instanton

[Potential]
potential = rgpot

[Instanton]
reactant_filename = reactant.con
product_filename = product.con
; start from a converged band instead of the straight line
initial_path = neb.con
beads = 256
springs = trotter
beta_hbar_omega = 30
force_tolerance = 1e-3
; one finite-difference Hessian every 4 beads, linear in between
hessian_stride = 4

springs = eco is refused. The instanton is a Trotter discretisation of the ring.

beta_hbar_omega sets the imaginary time in units of \(1/\omega\) of the stiffer minimum along the line between them. It must be long enough for the kink to relax into both wells. results.dat reports instanton_mode_separation, the second eigenvalue of \(J\) over its zero mode. Values above \(10^3\) mean the kink is isolated; small values mean \(\beta\hbar\) is too short. With no atom fixed, the product is aligned to the reactant first: its mass-weighted mean displacement is removed, and for a cluster its best rotation as well.

With that reactant rotation removed, each iteration evaluates every interior bead of the kink in one call. Under [RgpotPot] ranks_per_image, that call spreads the beads over the Car-Parrinello molecular dynamics (CPMD) calculator groups the same way a NEB spreads its images.

instanton.con writes one frame per bead, with imaginary_time_fs. The keys are:

Key	Meaning
tunnel_splitting_instanton	Splitting: \(\Delta_0\), eV
instanton_action	Action: \((S - S_\text{well})/\hbar\)
tls_energy_instanton	Energy: \(\sqrt{\Delta^2 + \Delta_0^2}\), eV
tunnel_asymmetry	Asymmetry: \(V(\text{product}) - V(\text{reactant})\), eV
instanton_temperature_K	Temperature: \(1/(k_B \beta)\) for the imaginary time used
instanton_mode_separation	Separation: how well the kink’s translation separates from the other modes
instanton_symmetric	Symmetry: 1 when {math}`\beta
instanton_beta_asymmetry	Magnitude: {math}`\beta

The propagator ratio measures the splitting \(\Delta_0\) when the two wells lie within a small fraction of \(k_B T\) of each other. instanton_symmetric = 0 flags a pair outside that window. The job still writes the path and the action, but no tunnel_splitting_instanton, and reports success: the flag says why. For such a pair, set mode = rate, give the saddle and a temperature below the crossover, and read the rate in the next sections.

Which path object

NEB images and ring-polymer beads are different objects. An image is a point on a path in configuration space between two minima. Its springs are fictitious. The parallel force is removed. A bead is one imaginary-time slice of a single quantum system. Its springs are physical, with stiffness fixed by the temperature and the number of beads.

The columns are:

Object	Points	Springs	What it returns
mode = splitting	open string between two minima, started from a band when one is present	Euclidean action, no tangent projection	tunnelling splitting when the wells are close in energy
mode = rate	closed ring through one saddle	same action, stiffness set by \(T\) and \(N\)	thermal rate: the ring below the crossover, the parabolic factor above it
Centroid potential of mean force (PMF)	one ring per image, centroid held on the image	sampled, not minimised	quantum free-energy barrier along the path
Harmonic centroid string	a ring at each image, optimised	local harmonic quantum correction	a free-energy estimate as good as that harmonic well

The first two rows are job = instanton. The centroid potential of mean force is a constrained path-integral molecular dynamics sample, one thermostatted ring per image. That sample is not an optimisation, and it does not belong in this job. A string of harmonically corrected rings is a different calculation again, and this page does not implement it.

The rate below the crossover

mode = rate reads the reactant and saddle_filename (default saddle.con). The instanton is a closed ring, a first-order saddle of the ring-polymer potential, with one negative eigenvalue and one zero eigenvalue that cycles the beads. Richardson and Althorpe [INST_RA09] give

\[k Z_r = \frac{1}{\beta_N \hbar} \sqrt{\frac{B_N}{2\pi \beta_N \hbar^2}} \prod_k' \frac{1}{\beta_N \hbar |\omega_k|} \exp(-\beta_N U_N).\]

\(B_N\) is the sum of squared steps around the ring, \(\omega_k^2\) are the eigenvalues of the mass-weighted ring Hessian, and the prime leaves out the cyclic zero mode and the rigid translations and rotations. \(Z_r\) is the harmonic ring-polymer partition function of the reactant. \(U_N\) is the ring-polymer potential, the bead potentials plus the springs.

The crossover temperature is \(T_c = \hbar \omega_b / (2\pi k_B)\), with \(\omega_b\) the imaginary frequency at the saddle. At or above \(T_c\) the ring collapses onto the saddle. Above \(T_c\) the rate is the parabolic barrier factor

\[\kappa = \frac{\pi T_c / T}{\sin(\pi T_c / T)}\]

times the classical harmonic transition-state theory rate from the same Hessians. \(\kappa\) tends to 1 at high temperature, which is the one-bead limit. The factor diverges as \(T\) approaches \(T_c\) from above, and at \(T_c\) the job stops. A large factor means the temperature is close to the crossover.

[Main]
job = instanton

[Instanton]
mode = rate
reactant_filename = reactant.con
saddle_filename = saddle.con
temperature = 5
beads = 256
hessian_stride = 1

temperature is in kelvin. It must be positive. Below \(T_c\) the job optimises the ring. Above \(T_c\) it writes the parabolic rate and does not optimise a ring. The default 0 means the temperature was not set, and mode = rate then refuses to run. beads defaults to 256, the same default as the splitting. A ring of 32 beads at 5 K does not resolve a stiff bond: the path integral starts to converge once the bead count exceeds \(\beta \hbar \omega\) of the stiffest mode. hessian_stride of 1 takes a Hessian on every bead. A larger stride keeps a Hessian on every stride-th bead and interpolates linearly between those anchors. That interpolation is an approximation: the rate formula uses the Hessian of every bead.

half_ring defaults to enabled. On an even bead count the potential is evaluated from one turning point to the other and copied onto the mirror. An odd count keeps every bead. energy_shift (default 0, in eV) is subtracted from every bead potential and from the reactant and saddle energies in the rate.

The search is an index-1 Newton step on the ring Hessian. Without an initial_path, the rate job traces a steepest-descent path out of the saddle along both signs of its unstable mode and seeds the ring from it by the period condition below. On LJ13 that path costs 586 gradient calls. Cooling a cosine ring from 0.85 of the crossover is the fallback when no path can be built; it finds the ring only where the ring grows continuously out of the saddle as the temperature drops. Where it does not, the search walks to a neighbouring saddle, so a converged ring with no bead on either side of the saddle’s dividing plane (the plane normal to its unstable mode) is refused and no rate is written. The step climbs the mode that overlaps the last climb and turns every other negative curvature downhill; the imaginary-time cycle and the rigid motions of the whole ring, rebuilt from the current beads, are held in place and left out of the step. A converged gradient is classified with finite-difference bead Hessians, since the Bofill blocks can carry negative curvatures the surface does not have; a second negative curvature that survives is a higher-index stationary ring, and the search steps down that mode. The ring Hessian is block cyclic tridiagonal in the beads, and every solve, determinant and inertia count goes through a block LU of the open chain plus a low-rank Woodbury correction for the closure, the cycle, the rigid modes and each eigenvector-following flip: \(O(N f^3)\) for \(f\) degrees of freedom, and the \(Nf \times Nf\) matrix is never formed, so the same step serves a seven-atom cluster and a 254-atom cell. The lowest ring modes come from Lanczos on matrix-vector products. The bead curvature blocks start from the saddle Hessian (initial_hessians = saddle, no force calls) and follow accepted moves with a Bofill update, rebuilt from finite differences up to three times when the trust radius reaches its floor; initial_hessians = finite_difference takes \(2f\) gradient calls per bead first. The rate uses the bead Hessians chosen by hessian_final, not that update. On the one-dimensional Eckart barrier the search converges in 4 to 7 steps from either seed.

temperatures is a comma-separated list in kelvin. The search starts at the highest and each ring starts the next, colder one. An empty list uses temperature. bead_ladder (default off) starts a ring of at least 16 beads at a quarter of that count and doubles. hessian_final is recomputed: the rate takes a Hessian on every stride-th bead.

On a rate calculation, initial_path is a band over the barrier. The ring starts on the closed orbit of that band whose period is \(\beta \hbar\), and the same band carries a one-dimensional WKB rate. rate_instanton.dat has one row per temperature, with columns T_K, T_c_K, beads, converged, iterations, U_N_eV, negative_modes, ln_k_per_s, k_per_s, ln_k_htst_per_s, barrier_effective_eV, ln_k_wkb_path_per_s, ln_k_parabolic_per_s and parabolic_factor. Above \(T_c\) the instanton columns are empty and the last two hold the parabolic rate. Below \(T_c\) those two are empty. A run at more than one temperature also writes instanton_<T>K.con. The coldest temperature is written to instanton.con.

With no atom fixed, both the reactant and the instanton omit the three translations. A rotation is omitted when it is a zero mode of the reactant Hessian, which a free cluster has and a crystal does not. A cluster in a large periodic cell is told apart by that Hessian, not by the periodic flag. The springs along those directions stay, so they cancel between the instanton and the reactant.

results.dat reports the rate. The keys are:

Key	Meaning
rate_instanton	\(k\) in s^{-1}
rate_instanton_log	\(\ln(k\,/\,\mathrm{s}^{-1})\)
rate_htst_log	classical harmonic transition-state theory, the same logarithm
parabolic_factor	\(\kappa\), above \(T_c\)
rate_parabolic	\(\kappa\) times the harmonic TST rate, s^{-1}
rate_parabolic_log	\(\ln(k\,/\,\mathrm{s}^{-1})\) of that rate
instanton_crossover_K	\(T_c\), K
instanton_negative_modes	negative eigenvalues of the ring Hessian; a first-order saddle has 1
instanton_zero_mode	the eigenvalue left out

Habershon et al. [INST_HMMM13] expect the sampled ring-polymer rate to lie within about a factor of two of the exact quantum rate between \(T_c\) and \(T_c/2\). That bound compares the sampled rate with the exact rate. It is not a comparison of this instanton with a free-energy profile.

The cubic metastable well, \(V = \omega_0^2 q^2/2 - g q^3/3\), is the check. Deep below the crossover its rate approaches the zero-temperature decay of Caldeira and Leggett [INST_CL83].

Path-integral quantum TST on planes

With pi_planes above 0, mode = rate follows the instanton with path-integral quantum transition-state theory (PI-QTST, Voth et al. [INST_VCM89]; review in Voth [INST_Vot93]). A ring polymer is sampled with its centroid held on each of a set of parallel planes. The centroid potential of mean force along the plane coordinate gives a free-energy barrier and a rate at every temperature in temperature or temperatures.

[Instanton]
mode = rate
saddle_filename = saddle.con
temperatures = 150, 105
pi_planes = 21
pi_beads = 32
pi_equilibration_steps = 500
pi_sampling_steps = 4000
pi_time_step = 0.5
pi_thermostat = pile
pi_direction = mode
pi_reactant_extent = 0.5

The coordinate. With \(q\) the mass-weighted displacement from the reactant over the free atoms and \(\hat n\) a unit vector in those coordinates, the plane coordinate is \(s = \hat n \cdot q\), in amu^0.5 Å. The reactant sits at \(s = 0\) and the saddle at \(s^* = \hat n \cdot q_\mathrm{saddle}\). pi_direction = mode (the default) takes \(\hat n\) from the unstable eigenvector of the saddle’s mass-weighted Hessian, oriented toward the saddle, so the last plane is the dividing surface normal to the barrier mode. When that mode makes more than 60 degrees with the reactant-saddle line, or the saddle has no negative eigenvalue, the job warns and uses pi_direction = line, the straight mass-weighted line from the reactant to the saddle. One normal serves every plane, so \(s\) is a linear coordinate and its mean force integrates to its free energy with no metric correction. The planes are fixed in the reactant’s frame. Translations of a structure with no atom fixed lie within every plane, because the unstable mode of the projected Hessian has no rigid component, but a rotation of a free cluster changes \(s\); fix an atom, or use a cell, when the sampling is long enough for the cluster to turn. The job warns when no atom is fixed in an aperiodic cell. pi_planes planes are spaced uniformly from \(s_0 = -\,\mathtt{pi\_reactant\_extent}\; s^*\), behind the reactant, to \(s^*\).

The sampling. On each plane the ring’s centroid starts where initial_path crosses the plane, or on the reactant-saddle line when no band is given. Projecting the centroid position and momentum holds it on the plane. The ring is carried from plane to plane: its centroid moves to the next plane and its internal modes keep their thermalised state. The free-ring springs are propagated exactly, so pi_time_step is limited by the physical vibrations, as for classical dynamics. pi_equilibration_steps are discarded, then pi_sampling_steps steps of pi_time_step fs record \(n \cdot f_c\), the centroid force along the plane normal. The bead forces of a step are one batch, so a calculator group carries the beads. pi_thermostat = pile puts PILE on the internal modes and a Langevin thermostat of time pi_pile_tau fs on the centroid within the plane, and pi_pile_scale (default 1) scales the critical damping of the internal modes. Below the crossover the lowest ring modes are soft at the barrier and overdamped at critical damping; on the Eckart check 0.5 cuts the mean-force error at 0.7 \(T_c\) by a third. piglet reads a normal-mode GLE from pi_gle_file, in the same format and with the same meaning as [Dynamics] path_gle_file. pi_seed seeds the noise.

The free energy. The mean force on the plane at \(s\) is

\[F'(s) = -\left\langle \hat n \cdot M^{-1/2} f_c \right\rangle_s ,\]

in eV per amu^0.5 Å, and \(F(s)\) is its trapezoid integral from \(s_0\). The production run is cut into ten equal blocks. The standard error of the block means is the error of \(F'(s)\). The errors of \(F\) and of the rate follow by linear propagation, with the planes taken as independent. The quantum free-energy barrier is \(\Delta F = F(s^*) - \min_{s < s^*} F(s)\). The log and results.dat give it beside the classical barrier \(V(\mathrm{saddle}) - V(\mathrm{reactant})\) and the instanton’s effective barrier.

The rate. With \(s\) a unit-mass coordinate,

\[k_\mathrm{PI\text{-}QTST} = \frac{1}{2}\sqrt{\frac{2}{\pi\beta}}\; \frac{e^{-\beta F(s^*)}}{\int_{s_0}^{s^*} e^{-\beta F(s)}\,ds},\]

with \(\beta = 1/k_B T\) in eV^{-1}. The prefactor is \(\tfrac12\langle|\dot s|\rangle\), half the mean speed of a free unit mass, in amu^0.5 Å per time unit of \(\sqrt{\mathrm{amu}\,\mathrm{Å}^2/\mathrm{eV}}\) (10.18 fs). The integral, again the trapezoid rule over the planes, is the reactant’s centroid density along \(s\). The job reports \(k\) in s^{-1}. The integral stops at the first plane, so that plane must sit several \(k_B T\) above the reactant minimum of \(F\); the job warns below 5 \(k_B T\). It also warns when \(F'(s^*)\) differs from zero by more than three standard errors: the maximum of the centroid free energy then lies off the classical saddle’s plane.

The formula is classical TST on the centroid free-energy surface. With one bead, \(F\) is the classical free energy and the rate is classical TST along \(s\). Near and above the crossover it carries the tunnelling correction of a symmetric barrier. For strongly asymmetric barriers well below the crossover it fails [INST_Vot93], and the instanton gives the rate there.

Output. rate_piqtst.dat has one row per temperature and plane, with columns T_K, s_amu05A, dF_ds_eV_per_amu05A, dF_ds_error, F_eV, F_error_eV and spread_max_A. piqtst_planes.con holds one frame per plane at the coldest temperature, at the production-averaged centroid. Its readcon spreads section holds each atom’s root-mean-square bead displacement from the centroid along x, y and z in Å. The scalars carry \(s\), \(F'\) and \(F\) with their errors, the temperature and the bead count. A run at more than one temperature also writes piqtst_planes_<T>K.con. results.dat gains, for the coldest temperature:

Key	Meaning
barrier_piqtst	\(\Delta F\), eV
barrier_piqtst_error	its standard error, eV
rate_piqtst	\(k_\mathrm{PI\text{-}QTST}\), s^{-1}
rate_piqtst_log	\(\ln(k\,/\,\mathrm{s}^{-1})\)
rate_piqtst_log_error	standard error of that logarithm
piqtst_s_star	\(s^*\), amu^0.5 Å
piqtst_dF_ds_star, piqtst_dF_ds_star_error	\(F'(s^*)\) and its error
piqtst_first_plane_kT	\(\beta (F(s_0) - \min F)\)
piqtst_temperature_K, piqtst_planes, piqtst_beads	the run

Cost: pi_planes times (pi_equilibration_steps plus pi_sampling_steps) steps per temperature, each step two force batches of pi_beads beads.

Recrossing and the RPMD rate

PI-QTST counts every centroid that reaches \(s^*\) moving forward as reactive. Some of those trajectories turn back. With pi_recrossing_parents above 0 the job measures that fraction, the Bennett-Chandler transmission factor of ring-polymer molecular dynamics (RPMD, Craig and Manolopoulos [INST_CM05]; the two-step scheme of Suleimanov et al. [INST_SAG13]), and reports

\[k_\mathrm{RPMD} = \kappa\, k_\mathrm{PI\text{-}QTST} .\]

[Instanton]
pi_recrossing_parents = 100
pi_recrossing_children = 20
pi_recrossing_time = 100
pi_recrossing_spacing = 50

The parents. A ring with its centroid held on the top plane \(s^*\), the dividing surface of the scan, is thermostatted as in the scan for pi_equilibration_steps. Every pi_recrossing_spacing steps after that its beads are one parent, pi_recrossing_parents in all.

The children. Each parent launches pi_recrossing_children momentum draws. A draw takes every ring normal mode from the Maxwell-Boltzmann distribution at \(\beta_P = \beta / P\) with the plane constraint removed, and runs twice, with \(p\) and with \(-p\), which halves the variance. Each child runs pi_recrossing_time fs of thermostat-free RPMD in steps of pi_time_step: velocity Verlet on the physical forces with the free ring propagated exactly in normal modes. Along the way the job records the centroid coordinate \(s(t)\). The centroid velocity along the coordinate is \(\dot s = \hat n \cdot M^{1/2} v_c\), and

\[\kappa(t) = \frac{\langle \dot s(0)\, h(s(t) - s^*) \rangle} {\langle \dot s(0)\, h(\dot s(0)) \rangle} ,\]

with \(h\) the step function. At \(t = 0\) the side is that of \(\dot s(0)\), the limit \(t \to 0^+\), so \(\kappa(0) = 1\). \(\kappa(t)\) falls as children recross and levels off once they have committed to a side. The reported \(\kappa\) is the mean of \(\kappa(t)\) over the last quarter of pi_recrossing_time. Its standard error is the jackknife over parents. Lengthen pi_recrossing_time when kappa_piqtst.dat has not levelled off by the last quarter.

What \(\kappa\) means. \(\kappa\) lies between 0 and 1. A value of 1 means no trajectory that crosses \(s^*\) forward returns, and PI-QTST is the RPMD rate. A value below 1 means the plane is not the dynamical bottleneck: either it is tilted from the barrier’s own dividing surface, or a bath coupling turns trajectories back. One bead gives the classical transmission through the plane. On a harmonic saddle whose plane normal makes an angle \(\theta\) with the unstable mode it is \(\sqrt{\cos^2\theta - \sin^2\theta\,\omega_b^2/\omega_\perp^2}\). The RPMD rate is independent of the choice of \(s^*\) in the long-time limit. The PI-QTST rate is not.

Output. kappa_piqtst.dat has the coldest temperature’s curve, columns t_fs and kappa. A run at more than one temperature also writes kappa_piqtst_<T>K.dat. rate_piqtst.dat gains the columns kappa, kappa_error, ln_k_rpmd_s and ln_k_rpmd_s_error, the same on every plane of a temperature. results.dat gains:

Key	Meaning
piqtst_kappa, piqtst_kappa_error	\(\kappa\) and its standard error
rate_rpmd	\(k_\mathrm{RPMD}\), s^{-1}
rate_rpmd_log	\(\ln(k_\mathrm{RPMD}\,/\,\mathrm{s}^{-1})\)
rate_rpmd_log_error	its standard error, both errors in quadrature

Option	Default	Meaning
pi_recrossing_parents	0	parent configurations; 0 is off, otherwise at least 2
pi_recrossing_children	20	momentum draws per parent, each run forward and reversed
pi_recrossing_time	100	fs per child, at least four pi_time_step
pi_recrossing_spacing	50	thermostatted steps between parents

Cost: pi_equilibration_steps plus pi_recrossing_parents times pi_recrossing_spacing constrained steps (two force batches each), and 2 * pi_recrossing_parents * pi_recrossing_children * pi_recrossing_time / pi_time_step child steps (one force batch each) per temperature.

Centroid and spread

Both modes also write instanton_centroid.con (and instanton_centroid_<T>K.con per temperature): one frame at the mean of the beads, with the readcon spreads section holding each atom’s root-mean-square displacement from that mean along x, y and z in Å. This is the delocalised configuration as a centroid plus a per-atom spread, the same representation a path-integral trajectory collapses to, so the atoms that tunnel and how far they spread read off one frame. Below the crossover the density is bimodal along the reaction path and the spread there is a width, not a Gaussian; the beads in instanton.con keep the full path. The frame carries spread_max, the largest entry, beside the temperature.

Checks

These checks use two Catch2 cases. The curved-valley case is Instanton splitting in a curved valley matches the exact gap. The corner case is The instanton cuts the corner the minimum energy path takes. The potential is \(V = V_0 (x^2 - 1)^2 + \tfrac{K}{2} (y - C (1 - x^2))^2\) at unit mass with \(K = 4\) eV/Ų. The exact gap comes from a fourth-order finite-difference Hamiltonian, converged to \(10^{-5}\).

\(V_0\) / eV	\(C\) / Å	exact \(\Delta_0\) / eV	instanton / exact	WKB on the valley floor / exact
0.12	0.35	2.520e-5	1.17	0.35
0.30	0.35	9.818e-8	1.08
0.12	0	2.040e-5	1.12	1.02
0.30	0	1.195e-7	1.07

The instanton’s error falls as the barrier deepens, the regime glass TLS sit in. The same cases tie the C++ path and splitting to an independent implementation of the discretisation to \(2 \times 10^{-3}\).

For the rate, The Eckart rate instanton matches the exact flux to its semiclassical error compares \(k Z_r\) through the symmetric Eckart barrier \(V_0 / \cosh^2(x/a)\) (\(V_0 = 0.425\) eV, \(a = 0.734\) amu^0.5 Å, \(T_c = 150\) K) with the exact flux \((2\pi\hbar)^{-1}\int P(E) e^{-\beta E} dE\) from Eckart’s transmission probability, at \(T = 0.5\,T_c\) and \(0.35\,T_c\):

beads	instanton / exact
64	0.94 to 0.96
128	0.93 to 0.94
\(N \to \infty\) (1/N² extrapolation)	0.928

In one dimension the instanton is the steepest-descent evaluation of the WKB thermal integral, so its limit shares the uniform WKB error; the Kemble integral along the path gives the same 0.928. The ring spectrum from the block chain matches the dense Hessian ties the chain’s determinant and inertia to a dense eigendecomposition, and A rigid mode leaves the instanton rate unchanged checks the rigid-mode bookkeeping through the search.

References

[INST_CL83]

A. Caldeira and A. Leggett. Quantum tunnelling in a dissipative system. Annals of Physics, 1983. doi:10.1016/0003-4916(83)90202-6.

[INST_CM05]

Ian R. Craig and David E. Manolopoulos. Chemical reaction rates from ring polymer molecular dynamics. The Journal of Chemical Physics, 122(8):084106, 2005. doi:10.1063/1.1850093.

[INST_HMMM13]

Scott Habershon, David E. Manolopoulos, Thomas E. Markland, and Thomas F. Miller. Ring-polymer molecular dynamics: quantum effects in chemical dynamics from classical trajectories in an extended phase space. Annual Review of Physical Chemistry, 64(1):387–413, April 2013. doi:10.1146/annurev-physchem-040412-110122.

[INST_RA09]

Jeremy O. Richardson and Stuart C. Althorpe. Ring-polymer molecular dynamics rate-theory in the deep-tunneling regime: connection with semiclassical instanton theory. The Journal of Chemical Physics, 131(21):214106, 2009. doi:10.1063/1.3267318.

[INST_SAG13]

Yury V. Suleimanov, Joshua W. Allen, and William H. Green. RPMDrate: bimolecular chemical reaction rates from ring polymer molecular dynamics. Computer Physics Communications, 184(3):833–840, 2013. doi:10.1016/j.cpc.2012.10.017.

[INST_Vot93] (1,2)

Gregory A. Voth. Feynman path integral formulation of quantum mechanical transition-state theory. The Journal of Physical Chemistry, 97(32):8365–8377, 1993. doi:10.1021/j100134a002.

[INST_VCM89]

Gregory A. Voth, David Chandler, and William H. Miller. Rigorous formulation of quantum transition state theory and its dynamical corrections. The Journal of Chemical Physics, 91(12):7749–7760, 1989. doi:10.1063/1.457242.