# Monte Carlo Integration

Monte Carlo integration relies on sampling from a certain distribution of random numbers. Importance sampling is a category of Monte Carlo approach which is widely used to compute multi-dimensional integral. In this approach, one first identifies a sampling distribution which resembles the target integrand. One then proceeds to sample the numbers from this distribution, and evaluates the modified integrand at these numbers. It is connected to umbrella sampling, and is popular among computational physicists. In importance sampling, the following approximation is made

$\int f(\mathbf{r}) g(\mathbf{r})\mathrm{d}\mathbf{r} \approx \frac{1}{L} \sum_{i=1}^{L}f(\mathbf{r}^i)$,

where $g(\mathbf{r})$ is a normalized sampling distribution. We use Monte Carlo technique to evaluate the classical flux-side correlation function, $C_{fs}(t)$. $C_{fs}(t)$ is used to compute classical rates, especially those that proceed at high temperatures or involve heavy reacting particles. $C_{\mathrm{fs}}(t)$ is defined as

$C_{\mathrm{fs}}(t) = \frac{1}{(2 \pi \hbar)^2}\int \mathrm{d}\mathbf{q}_0 \mathrm{d}\mathbf{p}_0\mathrm{e}^{-\beta H(\mathbf{p}_0,\mathbf{q}_0)}F(p_s,s)\Theta[s(\mathbf{q}_t)-s^{\ddagger}]$

$=\frac{1}{(2 \pi \hbar)^2}\int \mathrm{d}x_0\mathrm{d}y_0\mathrm{d}p_{x_0}\mathrm{d}p_{y_0} \mathrm{e}^{-\beta H(\mathbf{p}_0,\mathbf{q}_0)} \delta(x_0-s^{\ddagger})\frac{p_{x_0}}{m} \Theta[x_t-s^{\ddagger}]$

$=\frac{1}{(2 \pi \hbar)^2}\int \mathrm{d}x_0\mathrm{d}y_0\mathrm{d}p_{x_0}\mathrm{d}p_{y_0} \mathrm{e}^{- \frac{\beta(p_{x_0}^2+p_{y_0}^2)}{2 m}} \mathrm{e}^{-\beta V(x_0,y_0)}\delta(x_0-s^{\ddagger})\frac{p_{x_0}}{m} \Theta[x_t-s^{\ddagger}]$

where potential energy $V(x_0,y_0)=\frac{V_0}{\cosh^2{(x_0/a)}}+ \frac{1}{2} m \omega^2 y_0^2+ \frac{1}{2} m \omega^2 \sigma \mathrm{e}^{-\lambda(x_0-x_{\mathrm{shift}})^2 } y_0^2$, and $x_{\mathrm{shift}}$ represents the shift in coupling.

We want to evaluate the $C_{fs}(t)$ by using the method of Monte Carlo integration. For this, we first need to identify efficient sampling distribution. The normalized Gaussians from which we wish to sample $y_0$,$p_{x_0}$ and $p_{y_0}$ are defined as,

$P_1(y_0) = \omega\sqrt{\frac{\beta m }{2 \pi }} \mathrm{e}^{-\frac{1}{2} m \omega^2 y_0^2 \beta}$

$P_2(p_{x_0}) = \sqrt{\frac{\beta }{2 \pi m}} \mathrm{e}^{-\frac{1}{2}p_{x_0}^2 \frac{\beta}{m}}$

$P_3(p_{y_0}) = \sqrt{\frac{\beta }{2 \pi m}} \mathrm{e}^{-\frac{1}{2}p_{y_0}^2 \frac{\beta}{m}}$

We need not sample $x_0$, as this coordinate will be fixed by the properties of delta function. We therefore identify $g=P_1(y_0) P_2(p_{x_0}) P_3(p_{y_0})$. Therefore, after using the properties of delta function to integrate out $x_0$-degree of freedom, we can simply rewrite the integral as,

$C_{\mathrm{fs}}(t) \approx \frac{1}{L \beta \hbar^2 \omega} \sqrt{\frac{1}{ 2 \pi m \beta}} \mathrm{e}^{\frac{-\beta V_0}{\cosh^2{(s^{\ddagger}/a)}}} \sum_{i=1}^{L}p^{i}_{x_0} \mathrm{e}^{- \frac{\beta}{2} \sigma m \omega^2 (y_0^{i})^2 \mathrm{e}^{-(s^{\ddagger}-x_{\mathrm{shift}})^2 \lambda}} \Theta[x_t^{i}-s^{\ddagger}]$

$\Theta$ is the Heaviside step function, and $x_t$ is obtained by propagating the particle classically using velocity Verlet algorithm. The reaction coordinate $s$ is taken to be $x$-mode.

$C_{fs}(t)$ is computed for two different choices of dividing surface parameter $s^{\ddagger}$. In both cases, at long time, the value of $C_{fs}(t)/Z_r$ are identical. This suggests that classical rate is independent of the choice of dividing surface. The TST approximation to the classical rate is obtained in the short-time limit of $C_{fs}(t)$, which obviously depends on the choice of $s^{\ddagger}$.

Note: I prepared this problem for Advanced Kinetics exercise at ETH Zurich in Spring 2020.