Energy Decomposition along Reaction Coordinate and Energy Weighted Reactive Flux: I. Theory

Let $\vectorsym{x}(t)$ be an infinitely long trajectory of a system that consists of two stable states A and B, where $\vectorsym{x}$ represents the position coordinates of the system. In this long trajectory, the trajectory pieces that initiate from configurations in the state A and commit to the state B without returning to the state A are known to be reactive trajectories from A to B, which form the so-called TPE. We denote $U$ as an arbitrary ensemble, which could be the TPE or other non-equilibrium ensembles taken from the long trajectory and $\xi(\vectorsym{x})$ as a generalized coordinate. From a theoretical framework proposed in Ref. [²⁸] (here the notation and nomenclature in this reference are largely followed), the flux at a fixed location $\xi^{\prime}$ of $\xi(\vectorsym{x})$ in the ensemble $U$ is defined as

$$F^{U}(\xi^{\prime})=\int_{0}^{T}\dot{\xi}(\vectorsym{x}(t))\delta(\xi(\vectorsym{x}(t))-\xi^{\prime})h_{U}(t)dt$$

(1)

where $h_{U}(t)$ is an indicator function for the ensemble $U$, that is $h_{U}(t)=1$ if $\vectorsym{x}(t)$ belongs to $U$ and $h_{U}(t)=0$ otherwise.

If $\vectorsym{y}\equiv\{y_{1},y_{2},\cdots,y_{n}\}$ is a set of CVs, which are functions of $\vectorsym{x}(t)$, the current (from now on we use current when referring to a vector field and flux when a scalar is referred) in the space of CVs in the ensemble $U$ can then be defined as

$$\vectorsym{J}^{U}(\vectorsym{y})=\int_{0}^{T}\dot{\vectorsym{y}}(t)\delta(\vectorsym{y}(t)-\vectorsym{y})h_{U}(t)dt$$

(2)

Given $f(\vectorsym{x})$ being an arbitrary function of $\vectorsym{x}$, we consider a particular ensemble average of the quantity $f(\vectorsym{x})$ over the surface $\xi(\vectorsym{x})=\xi^{\prime}$ in the ensemble $U$,

$$\langle f(\vectorsym{x})\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}\equiv\frac{\int_{0}^{T}f(\vectorsym{x}(t))\dot{\xi}(\vectorsym{x}(t))\delta(\xi(\vectorsym{x}(t))-\xi^{\prime})h_{U}(t)dt}{\int_{0}^{T}\dot{\xi}(\vectorsym{x}(t))\delta(\xi(\vectorsym{x}(t))-\xi^{\prime})h_{U}(t)dt}$$

(3)

It is inspired by the principal curve proposed in the transition path theory ³⁵ and is actually the average over the surface $\xi(\vectorsym{x})=\xi^{\prime}$ of $f(\vectorsym{x})$ that is weighted with the current intensity passing through the point $\vectorsym{x}$ on the surface $\xi(\vectorsym{x})=\xi^{\prime}$, that is,

$$\langle f(\vectorsym{x})\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}=\frac{\int f(\vectorsym{x})\vectorsym{J}^{U}(\vectorsym{x})\cdot\nabla\xi(\vectorsym{x})\delta(\xi(\vectorsym{x})-\xi^{\prime})d\vectorsym{x}}{\int\vectorsym{J}^{U}(\vectorsym{x})\cdot\nabla\xi(\vectorsym{x})\delta(\xi(\vectorsym{x})-\xi^{\prime})d\vectorsym{x}}$$

(4)

Since a point on the principal curve given in Eq. (63) in Ref. [ ³⁵] can be viewed as a conditional canonical expectation of points in an isosurface along the principal curve itself under the assumptions: (1) the committor can be approximated by a function of the principal curve, (2) the curvature effects of the isosurface is negligible, and (3) the overdamped dynamics is assumed. In the following, we will show that $\langle f(\vectorsym{x})\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}$ is a conditional canonical expectation of $f(\vectorsym{x})$ under similar assumptions. For a diffusion process where the Smoluchowski equation is applicable, the reactive current is given by ³⁶

$$\vectorsym{J}_{R}(\vectorsym{x})=\rho(\vectorsym{x})\vectorsym{D}(\vectorsym{x})\nabla q(\vectorsym{x})$$

(5)

where $\rho(\vectorsym{x})$ is the equilibrium probability density in canonical ensemble, $\vectorsym{D}(\vectorsym{x})$ is the position-dependent diffusion matrix, and $q(\vectorsym{x})$ is the committor function. If $\xi(\vectorsym{x})$ is a monotonically increasing function of the committor, i.e., $q(\vectorsym{x})=g(\xi(\vectorsym{x}))$ and thus $\nabla q=g^{\prime}(\xi)\nabla\xi$, and $U$ is chosen to be the TPE (then $\vectorsym{J}^{U}(\vectorsym{x})$ is proportional to $\vectorsym{J}_{R}(\vectorsym{x})$), we have,

	$\displaystyle\langle f(\vectorsym{x})\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}=$	$\displaystyle\frac{\int f(\vectorsym{x})\rho(\vectorsym{x})g^{\prime}(\xi(\vectorsym{x}))\nabla\xi(\vectorsym{x})\cdot\vectorsym{D}(\vectorsym{x})\nabla\xi(\vectorsym{x})\delta(\xi(\vectorsym{x})-\xi^{\prime})d\vectorsym{x}}{\int\rho(\vectorsym{x})g^{\prime}(\xi(\vectorsym{x}))\nabla\xi(\vectorsym{x})\cdot\vectorsym{D}(\vectorsym{x})\nabla\xi(\vectorsym{x})\delta(\xi(\vectorsym{x})-\xi^{\prime})d\vectorsym{x}}$		(6a)
	$\displaystyle=$	$\displaystyle\frac{\int f(\vectorsym{x})\rho(\vectorsym{x})\nabla\xi(\vectorsym{x})\cdot\vectorsym{D}(\vectorsym{x})\nabla\xi(\vectorsym{x})\delta(\xi(\vectorsym{x})-\xi^{\prime})d\vectorsym{x}}{\int\rho(\vectorsym{x})\nabla\xi(\vectorsym{x})\cdot\vectorsym{D}(\vectorsym{x})\nabla\xi(\vectorsym{x})\delta(\xi(\vectorsym{x})-\xi^{\prime})d\vectorsym{x}}$		(6b)

Under the assumption of a local transition tube, in which most of the reactive current are located, the curvature of the surface $\xi(\vectorsym{x})=\xi^{\prime}$ in the transition tube could be negligible and the diffusion on it may be approximately position independent (not necessarily isotropic), then $\nabla\xi(\vectorsym{x})\cdot\vectorsym{D}(\vectorsym{x})\nabla\xi(\vectorsym{x})$ is approximately a constant and can be cancelled out in the above fraction. We thus arrive,

$$\langle f(\vectorsym{x})\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}\approx\frac{\int f(\vectorsym{x})\rho(\vectorsym{x})\delta(\xi(\vectorsym{x})-\xi^{\prime})d\vectorsym{x}}{\int\rho(\vectorsym{x})\delta(\xi(\vectorsym{x})-\xi^{\prime})d\vectorsym{x}}$$

(7)

Although we assume $f(\vectorsym{x})$ being a function of $\vectorsym{x}$, the weighted ensemble average can be generalized to any time series $f(t)$ extracted from trajectories and Eq. (3) is reformulated into

$$\langle f(t)\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}\equiv\frac{\int_{0}^{T}f(t)\dot{\xi}(\vectorsym{x}(t))\delta(\xi(\vectorsym{x}(t))-\xi^{\prime})h_{U}(t)dt}{\int_{0}^{T}\dot{\xi}(\vectorsym{x}(t))\delta(\xi(\vectorsym{x}(t))-\xi^{\prime})h_{U}(t)dt}$$

(8)

In case that $\xi(\vectorsym{x})$ is a good one-dimensional representation of the reaction coordinate and the flux along it ($F^{U}(\xi)$) presents a plateau region, in which $F^{U}(\xi)$ is a constant and equals the number of transition paths $N_{\rm T}$. We then have (Appendix A)

$$\frac{d\langle f(\vectorsym{x})\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}}{d\xi^{\prime}}=\langle\frac{df(\vectorsym{x}(t))}{dt}\frac{1}{\dot{\xi}(\vectorsym{x}(t))}\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}\quad\text{for $\xi^{\prime}$ in the plateau region}$$

(9)

In case that $\xi(\vectorsym{x})$ is the committor function or a monotone of the committor, Eq. (9) holds everywhere. We will show the usefulness of Eqs. (7) and (9) in the next subsection.

Given an independent and complete set of configurational coordinates $\vectorsym{z}\equiv\{z_{\rm 1},z_{\rm 2},\cdots,z_{\rm n}\}$, we consider a special case that $f(\vectorsym{z})$ is the potential energy $V(\vectorsym{z})$ and then

$$\langle\frac{dV(\vectorsym{z}(t))}{dt}\frac{1}{\dot{\xi}(t)}\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}=\sum_{i}\langle\frac{\partial V(\vectorsym{z}(t))}{\partial z_{\rm i}}\frac{\dot{z}_{\rm i}(t)}{\dot{\xi}(t)}\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}$$

(10)

Here, we avoid the explicit expressions of $\vectorsym{z}$ as functions of $\vectorsym{x}$ to simply the notations. Note that the dimension of $\frac{\partial V(\vectorsym{z}(t))}{\partial z_{\rm i}}\frac{\dot{z}_{\rm i}(t)}{\dot{\xi}(t)}$ is the same as the one of the generalized force corresponding to the generalized coordinate $\xi$, its integral over $\xi$ gives terms of the same dimension as energy, we thus formally defined the energy decomposed onto a coordinate $z_{\rm i}$ along the coordinate $\xi$ as

$$V_{z_{\rm i},\xi}=\int\langle\frac{\partial V(\vectorsym{z}(t))}{\partial z_{\rm i}}\frac{\dot{z}_{\rm i}(t)}{\dot{\xi}(t)}\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}d\xi^{\prime}$$

(11)

In the following, we will discuss the properties of the energy component in Eq. (11) in various cases under the conditions: (1) The ensemble $U$ is chosen to be the TPE; (2) The system can be described as a reaction coordinate in interaction with a large number of bath modes ³⁷; (3) The behaviours of bath modes in the TPE resemble the ones in the equilibrium ensemble.

The energy component in Eq. (11) can be related to the emergent potential energy proposed in Ref. [²⁵]. The emergent potential energy on a coordinate $z_{\rm i}$ along $\xi$ is defined as

$$\frac{\langle\Delta V_{z_{\rm i}}(\vectorsym{z};\xi^{\prime},\Delta\xi)\rangle}{\Delta\xi}=\langle\frac{\partial V(\vectorsym{z})}{\partial z_{\rm i}}\frac{\dot{z}_{\rm i}}{\dot{\xi}}\text{sgn}(\dot{\xi})\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}=\langle\frac{\partial V(\vectorsym{z})}{\partial z_{\rm i}}\frac{\dot{z}_{\rm i}}{|\dot{\xi}|}\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}$$

(20)

where sgn$(x)$ is the sign function and $\langle\cdots\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}$ stands for an ensemble average over the ensemble $U$ at $\xi(\vectorsym{z})=\xi^{\prime}$, which is defined as

$$\langle f(t)\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}\equiv\frac{\int_{0}^{T}f(t)\delta(\xi(\vectorsym{x}(t))-\xi^{\prime})h_{U}(t)dt}{\int_{0}^{T}\delta(\xi(\vectorsym{x}(t))-\xi^{\prime})h_{U}(t)dt}$$

(21)

On the other hand, the energy component on $z_{\rm i}$ along $\xi$ in this work can be rewritten with similar notations as

$$\langle\frac{\partial V(\vectorsym{z}(t))}{\partial z_{\rm i}}\frac{\dot{z}_{\rm i}(t)}{\dot{\xi}(t)}\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}=\frac{\langle\frac{\partial V(\vectorsym{z})}{\partial z_{\rm i}}\dot{z}_{\rm i}\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}}{\langle\dot{\xi}\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}}$$

(22)

Recall the definition of a transmission coefficient analogue suggested in Ref. [²⁸],

$$\check{\kappa}^{U}(\xi^{\prime})=\frac{\int_{0}^{T}\dot{\xi}(\vectorsym{x}(t))\delta(\xi(\vectorsym{x}(t))-\xi^{\prime})h_{U}(t)dt}{\int_{0}^{T}|\dot{\xi}(\vectorsym{x}(t))|\delta(\xi(\vectorsym{x}(t))-\xi^{\prime})h_{U}(t)dt}=\frac{\langle\dot{\xi}\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}}{\langle|\dot{\xi}|\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}}$$

(23)

The multiplication of the energy component in Eq. (22) and the transition coefficient analogue in Eq. (23) gives

$$\langle\frac{\partial V(\vectorsym{z}(t))}{\partial z_{\rm i}}\frac{\dot{z}_{\rm i}(t)}{\dot{\xi}(t)}\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}*\check{\kappa}^{U}(\xi^{\prime})=\frac{\langle\frac{\partial V(\vectorsym{z})}{\partial z_{\rm i}}\dot{z}_{\rm i}\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}}{\langle|\dot{\xi}|\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}}$$

(24)

When $\frac{\partial V(\vectorsym{z})}{\partial z_{\rm i}}\frac{\dot{z}_{\rm i}}{|\dot{\xi}|}$ and $|\dot{\xi}|$ are independent to each other in the ensemble $U$, that is,

$$\langle\frac{\partial V(\vectorsym{z})}{\partial z_{\rm i}}\dot{z}_{\rm i}\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}=\langle\frac{\partial V(\vectorsym{z})}{\partial z_{\rm i}}\frac{\dot{z}_{\rm i}}{|\dot{\xi}|}\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}\langle|\dot{\xi}|\rangle_{\xi(\vectorsym{z})=\xi^{\prime},U}$$

(25)

then the multiplication of the energy component and the transition coefficient analogue gives the emergent potential energy.

Although the emergent potential energy provide a rigorous way to quantify from an energetic viewpoint the relative importance of a coordinate for a rare transition when the TPE is available, the estimation of it suffers from numerical instabilities in the region where the probability density in the TPE is very low. Here, we proposed the following quantity to appraise the relevance of a coordinate to the reaction coordinate

		$\displaystyle\frac{d\hat{A}^{U}(\xi)}{d\xi}|_{\xi(\vectorsym{x})=\xi^{\prime}}\equiv\frac{d\tilde{A}^{U}(\xi)}{d\xi}|_{\xi(\vectorsym{x})=\xi^{\prime}}*\frac{F^{U}(\xi^{\prime})}{N_{\rm T}}=\frac{1}{N_{\rm T}}\int_{0}^{T}\frac{dV(\vectorsym{x}(t))}{dt}\delta(\xi(\vectorsym{x}(t))-\xi^{\prime})h_{U}(t)dt$		(26a)
	$\displaystyle\Rightarrow$	$\displaystyle\hat{A}^{U}(\xi_{\rm max})-\hat{A}^{U}(\xi_{\rm min})=\int^{\xi_{\rm max}}_{\xi_{\rm min}}\frac{d\tilde{A}^{U}(\xi)}{d\xi}|_{\xi(\vectorsym{x})=\xi^{\prime}}d\xi^{\prime}=\frac{1}{N_{\rm T}}\int_{0}^{T}\frac{dV(\vectorsym{x}(t))}{dt}h_{U}(t)dt$		(26b)

where, $\xi_{\rm min}$ and $\xi_{\rm max}$ are the lower and upper bounds of $\xi$ in the ensemble $U$, respectively. When the flux $F^{U}(\xi)$ along $\xi$ is vanishingly small, the term $\frac{F^{U}(\xi^{\prime})}{N_{\rm T}}$ will render $\frac{d\check{A}^{U}(\xi)}{d\xi}|_{\xi(\vectorsym{x})=\xi^{\prime}}$ to be close to zero and this quantity is unlikely to suffer from numerical instability. Since the difference of $\hat{A}^{U}(\xi)$ at the two boundaries of any coordinate are the same, it will be convenient via visual inspection to compare $\hat{A}^{U}(\xi)$ along different coordinates. When $\xi$ is the reaction coordinate, $\tilde{A}^{U}(\xi)$ is reproduced. When $\xi$ is a bath mode,

	$\displaystyle\frac{d\hat{A}^{U}(\xi)}{d\xi}|_{\xi(\vectorsym{x})=\xi^{\prime}}\approx\frac{T^{U}\rho^{U}(\xi^{\prime})}{N_{\rm T}}\langle\frac{dV(\vectorsym{x}(t))}{dq_{\rm s}}\dot{q}_{\rm s}\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}$		(27a)
	$\displaystyle\text{with}\quad T^{U}\equiv\int_{0}^{T}h_{U}(t)dt,\quad\rho^{U}(\xi^{\prime})\equiv\frac{\int_{0}^{T}\delta(\xi(\vectorsym{x}(t))-\xi^{\prime})h_{U}(t)dt}{\int_{0}^{T}h_{U}(t)dt}$		(27b)

where $T^{U}$ is the accumulated time in the ensemble $U$ and $\rho^{U}(\xi^{\prime})$ is the probability density at $\xi(\vectorsym{x})=\xi^{\prime}$ in the ensemble $U$. It should be reasonable to assume that $\langle\frac{dV(\vectorsym{x}(t))}{dq_{\rm s}}\dot{q}_{\rm s}\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}$ is a constant and then one immediately has

	$\displaystyle\hat{A}^{U}(\xi)\approx\overline{\Delta V^{U}}\int_{\xi_{\rm min}}^{\xi}\rho^{U}(\xi^{\prime})d\xi^{\prime}$		(28a)
	$\displaystyle\text{with}\quad\overline{\Delta V^{U}}\equiv\frac{1}{N_{\rm T}}\int_{0}^{T}\frac{dV(\vectorsym{x}(t))}{dt}h_{U}(t)dt$		(28b)

It is interesting to note that the average energy weighted with current intensity along the reaction coordinate is directly related to the free energy along it, as indicated in Eqs. (9) and (12). While the conditional canonical average of potential energy gives the internal energy, the entropic component of the free energy seems mainly arises from the curvature of the transition tube or from the shape change of transition tube intersecting with the isosurfaces of the reaction coordinate. It could be insightful in the configurational space to look at an energy weighted current, which is defined as

$$\vectorsym{J}_{V}(\vectorsym{x})\equiv\int_{0}^{T}V(\vectorsym{x}(t))\dot{\vectorsym{x}}(t)\delta(\vectorsym{x}(t)-\vectorsym{x})dt=V(\vectorsym{x})\vectorsym{J}(\vectorsym{x})$$

(29)

One immediately has,

$$\nabla\vectorsym{J}_{V}(\vectorsym{x})=\nabla V(\vectorsym{x})\cdot\vectorsym{J}(\vectorsym{x})+V(\vectorsym{x})\nabla\cdot\vectorsym{J}(\vectorsym{x})$$

(30)

In case of a steady state, i.e., $\nabla\cdot\vectorsym{J}(\vectorsym{x})=0$, then

$$\nabla\vectorsym{J}_{V}(\vectorsym{x})=\vectorsym{J}(\vectorsym{x})\cdot\nabla V(\vectorsym{x})=\int_{0}^{T}\frac{dV(\vectorsym{x}(t))}{dt}\delta(\vectorsym{x}(t)-\vectorsym{x})dt$$

(31)

Thus, the divergence of the energy weighted current equals the directional derivative of the potential energy in the direction of the reactive current for a steady state process. $\vectorsym{J}(\vectorsym{x})\cdot\nabla V(\vectorsym{x})$ can also be viewed as the flux intensity along the gradient of the potential energy and is weighted by the magnitude of the gradient.

The energy weighted current associated with an ensemble $U$ is then defined as,

$$\vectorsym{J}^{U}_{V}(\vectorsym{x})\equiv\int_{0}^{T}V(\vectorsym{x}(t))\dot{\vectorsym{x}}(t)\delta(\vectorsym{x}(t)-\vectorsym{x})h_{U}(t)dt$$

(32)

The projection of $\vectorsym{J}^{U}_{V}(\vectorsym{x})$ to the subspace of CVs gives the energy weighted current in the CVs subspace (Appendix B),

$$\vectorsym{J}^{U}_{V}(\vectorsym{y})=\int_{0}^{T}V(\vectorsym{x}(t))\dot{\vectorsym{y}}(t)\delta(\vectorsym{y}(t)-\vectorsym{y})h_{U}(t)dt$$

(33)

Inspired by Eq. (29), we thus seek for an energy function $\breve{A}^{U}_{V}(\vectorsym{y})$ of $\vectorsym{y}$ (thus a multidimensional free energy analogue in the subspace of CVs) such that $\vectorsym{J}^{U}_{V}(\vectorsym{y})=\breve{A}^{U}_{V}(\vectorsym{y})\vectorsym{J}^{U}(\vectorsym{y})$. One immediately finds that

$$\breve{A}^{U}_{V}(\vectorsym{y})=\frac{\vectorsym{J}^{U}_{V}(\vectorsym{y})\cdot\vectorsym{J}^{U}(\vectorsym{y})}{\vectorsym{J}^{U}(\vectorsym{y})\cdot\vectorsym{J}^{U}(\vectorsym{y})}$$

(34)

and that in the one-dimensional case $\breve{A}^{U}_{V}(\xi^{\prime})=\langle V(\vectorsym{x})\rangle_{\xi(\vectorsym{x})=\xi^{\prime},U}^{Flux}$.