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.

Cfs

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.