MP/SOFT Method

The Matching-Pursuit/Split Operator Fourier Transform (MP/SOFT) method is essentially the SOFT approach implemented in coherent-state representation, i.e., where the grid-based representation of $\tilde{\rho}({\bf x}_j,{\bf x'}_k;{\epsilon})$ is substituted by coherent-state expansions generated according to the Matching-Pursuit algorithm.

The MP/SOFT propagation of the initial condition entails the following steps:

<>Step [1]: Decompose $\tilde{\rho}({\bf x},{\bf x}';{\epsilon}) \equiv e^{-i V_0({\bf x}) \tau/2} \rho({\bf x},{\bf x}';{\epsilon})$ in a matching-pursuit coherent-state expansion:

$$\tilde{\rho}({\bf x},{\bf x}';{\epsilon}) \approx \sum_{j=1}^n c_{j} \phi_j( {\bf x}) [\phi'_j({\bf x}')]^*,$$ (540)


where $\phi_j({\bf x})$ and $\phi'_j({\bf x})$ are $N$-dimensional coherent-states defined as follows,

$$\phi_j({\bf x}) \equiv \prod_{k=1}^N A_{\phi_j}(k) e^{-\gamma_{\p... ...hspace{.1cm} e^{i \hspace{.1cm} p_{\phi_j}(k) \left(x(k)-x_{\phi_j}(k)\right)},$$ (541)


with complex-valued coordinates $x_{\phi_j}(k) \equiv r_{\phi_j}(k)+i d_{\phi_j}(k)$, momenta $p_{\phi_j}(k) \equiv g_{\phi_j}(k)+i f_{\phi_j}(k)$ and scaling parameters $\gamma_{\phi_j}(k) \equiv a_{\phi_j}(k)+i b_{\phi_j}(k)$. The normalization constants are $A_{\phi_j}(k) \equiv (a_{\phi_j}(k)/\pi)^{1/4}$   exp$[- \frac{1}{2} a_{\phi_j}(k) d_{\phi_j}(k)^2 - d_{\phi_j}(k) g_{\phi_j}(k) - \left(b_{\phi_j}(k) d_{\phi_j}(k) + f_{\phi_j}(k) \right)^2/(2 a_{\phi_j}(k))]$.

The expansion coefficients, introduced by Eq. (540), are defined as follows:

$$\begin{split} c_j \equiv \begin{cases} I_j, & \text{when}... ...le, & \text{for} \hspace{.3cm} j=2-n, \end{cases} \end{split}$$ (542)


where the overlap integral $I_j$ is defined as follows,

$$I_j \equiv \int d {\bf x}' d {\bf x} \hspace{.1cm} \phi_j({\bf x}) \tilde{\rho}({\bf x},{\bf x}';{\epsilon}) [\phi'_j({\bf x}')]^*.$$ (543)



• Step [2]: Analytically Fourier transform the coherent-state expansion to the momentum representation, apply the kinetic energy part of the Trotter expansion and analytically inverse Fourier transform the resulting expression back to the coordinate representation to obtain the imaginary-time evolved Boltzmann-operator matrix elements:

  $$\rho({\bf x},{\bf x}';{\epsilon} + i \tau) = \sum_{j=1}^n c_j e^{-i V_0({\bf x}) \tau/2} \widetilde{\phi}_j({\bf x}) [\phi'_j({\bf x}')]^*,$$ (544)

  
  
  where

  $$\widetilde{\phi}_j({\bf x}) \equiv \prod_{k=1}^N A_{\widetilde{\phi}_j}(k) \sqrt{\frac{m}{m+i \tau \gamma_{\widetilde{\phi}_j}(k)}} \left( \frac{\left( \frac{p_{\widetilde{\phi}_j}(k)}{\gamma_{\wid... ...- \frac{p_{\widetilde{\phi}_j}(k)^2}{2 \gamma_{\widetilde{\phi}_j}(k)} \right).$$ (545)

  
  

Note that the MP/SOFT approach reduces the computational task necessary for the imaginary- or real-time propagation of the Boltzmann operator matrix elements $\rho({\bf x},{\bf x}';{\beta})$ to the problem of recursively generating the coherent-state expansions introduced by Eq. (540). Coherent-state expansions are obtained by combining the matching pursuit algorithm and a gradient-based optimization method as follows:

• Step [1.1]. Evolve the complex-valued parameters, that define the initial trial coherent-states $\phi_j({\bf x})$ and $\phi'_j({\bf x})$, to locally maximize the overlap integral $I_j$, introduced in Eq. (543). The parameters $x_{\phi_1}(k), p_{\phi_1}(k), \gamma_{\phi_1}(k)$ and $x_{\phi'_1}(k), p_{\phi'_1}(k), \gamma_{\phi'_1}(k)$ of the corresponding local maximum define the first pair of coherent-states $\phi_1$ and $\phi'_1$ in the expansion introduced by Eq. (540) and the first expansion coefficient $c_1$, as follows: $\tilde{\rho}({\bf x},{\bf x}';\epsilon) = c_1 \phi_1({\bf x}) [\phi'_1({\bf x}')]^* + \varepsilon_1({\bf x},{\bf x}'),$ where $c_1 \equiv I_1$, as defined according to Eq. (543). Note that due to the definition of $c_1$, the residue $\varepsilon_1({\bf x},{\bf x}')$ does not overlap with the product state $\phi_1({\bf x}) [\phi'_1({\bf x}')]^*$. Therefore, the norm of the remaining residue $\varepsilon_1({\bf x},{\bf x}')$ is smaller than the norm of the initial target state $\tilde{\rho}({\bf x},{\bf x}';\epsilon)$ --i.e., $\parallel \varepsilon_1 \parallel < \parallel \tilde{\rho} \parallel$.
• Step [1.2]. Goto [1.1], replacing $\tilde{\rho}({\bf x},{\bf x}';{\epsilon})$ by $\varepsilon_1({\bf x},{\bf x}')$ --i.e., sub-decompose the residue by its projection along the direction of its locally optimum match as follows: $\varepsilon_1({\bf x},{\bf x}') = c_2 \phi_2({\bf x}) [\phi'_2({\bf x}')]^* + \varepsilon_2({\bf x},{\bf x}'),$ where

  $$c_2 \equiv \int d {\bf x}' d {\bf x} \hspace{.1cm} \phi_2({\bf x}) \varepsilon_1({\bf x},{\bf x}') [\phi'_2({\bf x}')]^*.$$ (546)

  
  
  Note that $\parallel \varepsilon_2 \parallel < \parallel \varepsilon_1 \parallel$, since $\varepsilon_2({\bf x},{\bf x}')$ is orthogonal to the product state $\phi_2({\bf x}) [\phi'_2({\bf x}')]^*$.

Step [1.2] is repeated each time on the resulting residue. After $n$ successive projections, the norm of the residue $\varepsilon_n$ is smaller than a desired precision $\epsilon$ --i.e., $\parallel \varepsilon_n \parallel = (1 - \sum_{j=1}^n \vert c_j \vert ^2)^{1/2} < \epsilon,$ and the resulting expansion is given by Eq. (540).
Note that norm conservation of $\hat{\tilde{\rho}}_{\epsilon}$ is maintained within a desired precision, just as in a linear orthogonal decomposition, although the coherent-states in the expansion are non-orthogonal basis-functions. It is important to mention that the computational bottleneck of the MP/SOFT method involves the calculation of overlap matrix elements $\langle \phi_j \mid e^{-i V_j(\hat{{\bf x}}) \tau/2} \mid \widetilde{\phi}_k \rangle$ and $\langle \phi_j \mid e^{-i V_j(\hat{{\bf x}}) \tau/2} \mid \phi_k \rangle$, where $\vert \phi_k \rangle$ and $\vert \widetilde{\phi}_k \rangle$ are localized Gaussians introduced by Eqs. (541) and (545), respectively. The underlying computational task is however trivially parallelized.

The overlap integrals are most efficiently computed in applications to reaction surface Hamiltonians where a large number of harmonic modes can be arbitrarily coupled to a few reaction (tunneling) coordinates. For such systems, the Gaussian integrals over harmonic coordinates can be analytically computed and the remaining integrals over reaction coordinates are efficiently obtained according to numerical quadrature techniques. For more general Hamiltonians, the overlap matrix elements can be approximated by analytic Gaussian integrals when the choice of width parameters $\gamma_j(k)$ allows for a local expansion of $V_j(\hat{{\bf x}})$ to second order accuracy. Otherwise, the quadratic approximation is useful for numerically computing the corresponding full-dimensional integrals according to variance-reduction Monte Carlo techniques.