Molecular Flow Monte Carlo

TL;DR

This paper introduces a novel continuous normalizing flow method for Monte Carlo molecular simulations that efficiently generates independent, low-energy configurations by transforming distributions based on two-body interactions, improving sampling efficiency.

Contribution

It presents a new continuous normalizing flow approach tailored for molecular systems, enabling efficient sampling by avoiding high-energy states and improving proposal distributions.

Findings

The inverse square flow effectively transforms distributions for molecular sampling.

The method improves Metropolis-Hastings efficiency in simple molecular systems.

It avoids extremely high energy configurations during sampling.

Abstract

In this paper, we suggest a novel sampling method for Monte Carlo molecular simulations. In order to perform efficient sampling of molecular systems, it is advantageous to avoid extremely high energy configurations while also retaining the ability to quickly generate new and independent trial states. Thus, we introduce a continuous normalizing flow method which can quickly generate independent states for various proposal distributions using a first-order differential equation. We define this continuous normalizing molecular flow approach based on two-body intermolecular interactions to achieve a probability distribution transformation method which yields distributions which have probability densities of zero when molecule pairs are in close proximity; while in all other cases, the probability density is compressed such that it is spatial uniform. This transform provides the proposal distribution which generates no states of extremely high energy. We find that an inverse square flow is applicable as the continuous normalizing molecular flow. Using the transformed distribution, we can perform the Metropolis-Hastings method more efficiently. The high efficiency of the proposed method is demonstrated using simple molecular systems.

Equations

$$p({\mathbf{r}}^N)\propto \exp \{-\beta U({\mathbf{r}}^N)\},$$

$$\mathrm{where}U({\mathbf{r}}^N)=\sum _{i<j}{U}_{ij}(\Vert {\mathbf{r}}_i-{\mathbf{r}}_j\Vert ),$$

$$\frac{d\mathbf{z}}{dt}=f(\mathbf{z}(t),t),$$

$$\frac{\partial \log p(\mathbf{z}(t))}{\partial t}=-\mathrm{tr}\left(\frac{df}{d\mathbf{z}(t)}\right)$$

$$\displaystyle\frac{d{\bm{r}}_{i}}{dt}=\sum_{i\neq j}{f_{ij}(\|{\bm{r}}_{i}-{\bm{r}}_{j}\|)({\bm{r}}_{i}-{\bm{r}}_{j}})\$$

$$\mathrm{tr}\left(\frac{df}{d{\mathbf{r}}^N(t)}\right)=\sum _{i,j}D{f}_{ij}(r_{ij})+r_{ij}\frac{d{f}_{ij}}{d{r}_{ij}}(r_{ij}),$$

$$f_{ij}(r)=\frac{1}{r^D}.$$

$$\frac{d{\mathbf{r}}_i}{dt}=\sum _{\mathbf{n}}\sum _{i\ne j}{f}_{ij}(\Vert {\mathbf{r}}_i-{\mathbf{r}}_j-\mathbf{n}\Vert )({\mathbf{r}}_i-{\mathbf{r}}_j-\mathbf{n})$$

$$\mathrm{tr}\left(\frac{df}{d{\mathbf{r}}^N(t)}\right)=\sum _{\mathbf{n}}\sum _{i\ne j}D{f}_{ij}(r_{ij,\mathbf{n}})+r_{ij,\mathbf{n}}\frac{d{f}_{ij}}{d{r}_{ij,\mathbf{n}}}(r_{ij,\mathbf{n}}),$$

$$\sum _{\mathbf{n}}\Vert \mathbf{r}-\mathbf{n}{\Vert }^{-3}(\mathbf{r}-\mathbf{n})=\sum _{\mathbf{n}}{\nabla }_{\mathbf{r}}\frac{1}{\Vert \mathbf{r}-\mathbf{n}\Vert }=$$

$${\nabla }_{\mathbf{r}}\left(\sum _{\mathbf{n}}\frac{\mathrm{erfc}(G\Vert \mathbf{r}-\mathbf{n}\Vert )}{\Vert \mathbf{r}-\mathbf{n}\Vert }+\frac{1-\mathrm{erfc}(G\Vert \mathbf{r}-\mathbf{n}\Vert )}{\Vert \mathbf{r}-\mathbf{n}\Vert }\right).$$

$$\sum _{\mathbf{n}}\Vert \mathbf{r}-\mathbf{n}{\Vert }^{-3}(\mathbf{r}-\mathbf{n})=$$

$${\nabla }_{\mathbf{r}}\left(\sum _{\mathbf{n}}\frac{\mathrm{erfc}(G\Vert \mathbf{r}-\mathbf{n}\Vert )}{\Vert \mathbf{r}-\mathbf{n}\Vert }+4\pi \sum _{\mathbf{k}\ne 0}\frac{e^{-\Vert \mathbf{k}{\Vert }^2/4G^2}{e}^{i\mathbf{k}\cdot \mathbf{r}}}{\Vert \mathbf{k}{\Vert }^2}\right),$$

$$\frac{d{\mathbf{r}}_i}{dt}=\sum _{i\ne j}{\mathbf{f}}^p({\mathbf{r}}_i-{\mathbf{r}}_j).$$

$$A({\mathbf{r}}_{\ast }^N,{\mathbf{r}}^N)=\min \left(1,\frac{p({\mathbf{r}}_{\ast }^N)}{p({\mathbf{r}}^N)}\frac{q({\mathbf{r}}^N|{\mathbf{r}}_{\ast }^N)}{q({\mathbf{r}}_{\ast }^N|{\mathbf{r}}^N)}\right),$$

$$\displaystyle q({\bm{r}}_{i*}|{\bm{r}}_{i})=\frac{1}{(\sqrt{2\pi}{r_{move}})^{3}}\exp\biggl{\{}-\frac{\|{\bm{r}}_{i*}-{\bm{r}}_{i}\|^{2}}{2{r_{move}}^{2}}\biggr{\}},$$

$$g(r)=\frac{⟨n(r)⟩}{4\pi r^{2}dr\rho },$$

$$U_{ij}(r)=4\left(r^{-12}-r^{-6}\right),$$

$${\nabla }_{\mathbf{r}}\left(\sum _{\mathbf{n}}\frac{\mathrm{erfc}(G\Vert \mathbf{r}-\mathbf{n}\Vert )}{\Vert \mathbf{r}-\mathbf{n}\Vert }+4\pi \sum _{\mathbf{k}\ne 0}\frac{e^{-\Vert \mathbf{k}{\Vert }^2/4G^2}{e}^{i\mathbf{k}\cdot \mathbf{r}}}{\Vert \mathbf{k}{\Vert }^2}\right)$$

$$\displaystyle-\sum_{\bm{n}}{\frac{{\bm{r}}-{\bm{n}}}{\|{\bm{r}}-{\bm{n}}\|^{3}}}\biggl{\{}{\rm erfc}(G\|{\bm{r}}-{\bm{n}}\|)$$

$$+\frac{2G}{\sqrt{\pi }}\Vert \mathbf{r}-\mathbf{n}\Vert e^{-G^2\Vert \mathbf{r}-\mathbf{n}{\Vert }^2}\}$$

$$-4\pi \sum _{\mathbf{k}\ne 0}\frac{e^{-\Vert \mathbf{k}{\Vert }^2/4G^2}}{\Vert \mathbf{k}{\Vert }^2}\mathbf{k}\sin (\mathbf{k}\cdot \mathbf{r}),$$

$$\frac{d{\mathbf{r}}_i}{dt}=\sum _{i\ne j}{\mathbf{f}}^p({\mathbf{r}}_i-{\mathbf{r}}_j).$$

$${\nabla }_{\mathbf{r}}\left(\sum _{\mathbf{n}}\frac{\mathrm{erfc}(G\Vert \mathbf{r}-\mathbf{n}\Vert )}{\Vert \mathbf{r}-\mathbf{n}\Vert }\right).$$

$$\mathrm{tr}\left(\frac{df}{d{\mathbf{r}}^N(t)}\right)=\sum _{i\ne j}\frac{2}{r_{ij}}\frac{d{h}_{ij}}{d{r}_{ij}}+\frac{d^2{h}_{ij}}{{d{r}_{ij}}^2},$$

$$\mathrm{tr}\left(\frac{df}{d{\mathbf{r}}^N(t)}\right)=\sum _{i\ne j}\frac{4G^2}{\sqrt{\pi }}\sum _{\mathbf{n}}{e}^{-G^2\Vert {\mathbf{r}}_i-{\mathbf{r}}_j-\mathbf{n}{\Vert }^2}$$

$$h_{ij}(r)=\sum _{\mathbf{n}}\frac{\mathrm{erfc}(G\Vert \mathbf{r}-\mathbf{n}\Vert )}{\Vert \mathbf{r}-\mathbf{n}\Vert }.$$

$$\mathrm{tr}\left(\frac{df}{d{\mathbf{r}}^N(t)}\right)=\sum _{i\ne j}\underset{G\to 0}{\lim }\frac{4G^2}{\sqrt{\pi }}\sum _{\mathbf{n}}{e}^{-G^2\Vert {\mathbf{r}}_i-{\mathbf{r}}_j-\mathbf{n}{\Vert }^2}$$

$$=\sum _{i\ne j}\frac{4}{\sqrt{\pi }}{\int }_{\mathbf{x}\in {\mathbb{R}}^3}{e}^{-\Vert \mathbf{x}{\Vert }^2}$$

$$=\sum _{i\ne j}\frac{4}{\sqrt{\pi }}{\pi }^{3/2}=\sum _{i\ne j}4\pi$$

Full text

1. Introduction

Monte Carlo molecular simulation is one of the most commonly used and powerful tools with which to investigate the physical properties of molecular systems. It has seen widespread use in various fields, e.g., polymers (?), proteins (?), nucleation (?), and hydrates (?). The method works by generating states that follow the Boltzmann distribution for specific molecular systems, for example,

[TABLE]

${\bm{r}}^{N}=({\bm{r}}_{1},{\bm{r}}_{2},...,{\bm{r}}_{N})$ denotes the total set of 3-dimensional $N$ molecule coordinates, $\beta$ is the inverse temperature, $U({\bm{r}}^{N})$ is the total potential energy of system, and $U_{ij}$ is the two-body intermolecular potential energy between molecules $i$ and $j$. To sample the Boltzmann distribution, we usually use dedicated Monte Carlo methods such as the Metropolis-Hastings algorithm (?), replica exchange method (?), and umbrella sampling (?), according to the characteristics of the system being studied. Despite the applicability of the Monte Carlo approach, performing large-scale simulations is not practical because existing Monte Carlo methods generate new states which are quite similar to preceding ones. The implication of this is that, for large systems, it takes a long time to generate a sufficient number of independent states, which are required to obtain unbiased results. There is, therefore, a need for fast and efficient methods with which to generate independent states for large systems.

Recently, a notable sampling method, continuous normalizing flow (?), was proposed. This method quickly generates independent states which follow various distributions using a parametric first-order differential equation. Additionally, the probability density of each state may also be calculated. The method is suitable for use with Monte Carlo simulations and has great potential to improve sampling efficiency.

In this paper, we examine the applicability of continuous normalizing flow for molecular systems. We define continuous normalizing molecular flow based on two-body intermolecular interactions. Furthermore, we find that inverse square flow is suitable for continuous normalizing molecular flow. This flow is a space-constant instantaneous change flow in 3-dimensional space, which can be used even in the presence of periodic boundary conditions. Using this flow as a sampler for proposal distributions, we present a highly efficient Metropolis-Hastings method, named Molecular Flow Metropolis-Hastings (MFMH). The efficiency of MFMH is evaluated through simulations of simple molecular systems.

2. Continuous normalizing molecular flow

Continuous normalizing flow generates the state ${\bm{z}}(t)$ by solving the first-order differential equation

[TABLE]

where $f$ is a uniform Lipschitz continuous function in ${\bm{z}}$, while also being continuous in $t$. The distribution $p({\bm{z}}(t))$ that state ${\bm{z}}(t)$ follows is easily calculated by solving following differential equation

[TABLE]

from the initial distribution $p({\bm{z}}(0))$. The term ${\rm tr}(df/d{\bm{z}})$ is referred to as the instantaneous change of variables. Here, let us consider ${\bm{z}}$ as the total set of molecular coordinates ${\bm{r}}^{N}=({\bm{r}}_{1},{\bm{r}}_{2},...,{\bm{r}}_{N})$ and $f$ as two-body intermolecular “force”. Solving first-order the differential equation gives

[TABLE]

for all $i$, where $f_{ij}(r)$ is the two-body intermolecular function. In this case, the instantaneous change of ${\bm{r}}^{N}$ is described as

[TABLE]

where $r_{ij}=\|{\bm{r}}_{i}(t)-{\bm{r}}_{j}(t)\|$, and $D$ is the number of spatial dimensions of the system (usually $D=3$). Note that if we regard continuous normalizing molecular flow as machine learning architecture, then this architecture has novel properties: (particle) permutation invariance and translational/rotational invariance in the presence of periodic boundaries.

3. Inverse square flow

Inverse square flow on no boundary condition

In this section, we introduce the inverse square flow which achieves a space-constant instantaneous change of ${\bm{r}}^{N}$ while preventing the generation of states with molecule pairs in close proximity. By solving the equation $Df+rdf/dr=0$, to satisfy the condition that the instantaneous change of ${\bm{r}}^{N}$ (Eq.6) must always be zero, we obtain the following solution:

[TABLE]

We may call this solution $D$-dimensional inverse square flow because the term ${\bm{r}}/r^{D}$ takes the same form as the inverse-square law. In Fig.1 we present the numerical results of a 2-dimensional inverse square flow. We see the distribution of a molecule transformed from a uniform distribution in response to the inverse square force from immobile molecules placed at the four corners of the plotted area. In the transformed distribution, we can confirm that the occurrence of molecule pairs positioned in close proximity is prevented, while the probability density is the same as the untransformed distribution since the instantaneous change of ${\bm{r}}^{N}$ is 0. Therefore, by using inverse square flow, we can create a transformation that does not change the probability density and prevents the generation of states where molecules are in close proximity. In Monte Carlo molecular simulation, any state ${\bm{r}}^{N}$ which includes molecule pairs in close proximity will always have a very high energy $U({\bm{r}}^{N})$. It is computationally inefficient to attempt the sampling of such states. Thus, the method presented here has the potential to greatly enhance the efficiency of sampling within Monte Carlo simulations. This promises to be especially useful for systems with strongly repulsive potentials and/or high-density phases.

Inverse square flow with periodic boundary conditions

We perform molecular simulations in the presence of periodic boundary conditions. With periodic boundary conditions, each molecule is affected by other molecules in the original cell and the images of molecules in the neighbouring periodic cells that surround the simulated cell (see Fig.2). For simplicity, let us consider that a cell is a cube with each side being a unit in length. Continuous normalizing molecular flow for periodic boundary conditions is described as

[TABLE]

and the instantaneous change of ${\bm{r}}^{N}$ is described as

[TABLE]

where $r_{ij,{\bm{n}}}=\|{\bm{r}}_{i}-{\bm{r}}_{j}-{\bm{n}}\|$ and ${\bm{n}}\in\mathbb{Z}^{3}$. Unfortunately, when simply applying inverse square flow to systems with periodic boundary conditions, the sum of all periodic forces (the righthand term in Eq.8) does not converge because the inverse square force is a long-range force. To handle inverse square flow, we use Ewald summation (?) for the calculation of this force when $D=3$. In Ewald summation, the inverse square force is divided into

[TABLE]

The left term and right term in Eq.10 are summed in real space and reciprocal space, respectively. $G$ is the screening factor that determines the relative proportion of the real and reciprocal space sums. Using a pairwise form of the Ewald sum (?), instead of direct calculation of the periodic sum, gives

[TABLE]

where ${\bm{k}}\in{2\pi\mathbb{Z}^{3}}$ is the reciprocal vector. This leads to converge for any $G\in(0,\infty)$. A non-derivative form of the pairwise force is given in Appendix A. By expressing Eq.11 as ${\bm{f}}^{p}({\bm{r}})$, the inverse square flow for periodic boundary conditions may be written in a more compact form:

[TABLE]

Note that this equation does not have the problem of non-convergence of the periodic sum. Notably, though we only use the inverse square flow, the periodic sum of the inverse square flow has a non-zero instantaneous change. The instantaneous change of the periodic sum of the inverse square flow is also space-constant, but the value is $\sum_{i\neq j}4\pi$ (a proof of this is given in Appendix B).

4. Molecular Flow Monte Carlo

Standard Metropolis-Hastings method

Let us consider the case where $p({\bm{r}}^{N})$ is the target distribution, such as in Eq.1 when states from $p({\bm{r}}^{N})$ are required. Although the mathematical expression of $p({\bm{r}}^{N})$ is usually well-known, the direct generation of states from $p({\bm{r}}^{N})$ is typically difficult. The Metropolis-Hastings (MH) method generates states by stochastically accepting states from a Markov chain. If we iteratively generate states from the Markov chain $q({\bm{r}}^{N}_{*}|{\bm{r}}^{N})$ (the proposal distribution) and accept states with the following probability:

[TABLE]

then the generated states will follow $p({\bm{r}}^{N})$ (?). In the MH method, the mathematical expression of $q({\bm{r}}^{N}_{*}|{\bm{r}}^{N})$ is required for the calculation of the acceptance probability. Or, if the proposal distribution is symmetric ($q({\bm{r}}^{N}_{*}|{\bm{r}}^{N})=q({\bm{r}}^{N}|{\bm{r}}^{N}_{*})$), the acceptance probability $A({\bm{r}}^{N}_{*},{\bm{r}}^{N})$ reduces to $\min\left(1,p({\bm{r}}^{N}_{*})/p({\bm{r}}^{N})\right)$, thus, under these conditions, a mathematical expression for $q({\bm{r}}^{N}_{*}|{\bm{r}}^{N})$ is not actually required in practice.

In Monte Carlo molecular simulation, we often move only one molecule at a time. The movement of many molecules often results in a large energy change, and the acceptance probability becomes exponentially small. To maintain a sufficiently high acceptance ratio, we only move one randomly selected molecule at a time.

The most commonly used proposal distribution for $q({\bm{r}}_{i*}|{\bm{r}}_{i})$ is a normal distribution:

[TABLE]

where ${r_{move}}\in\mathbb{R}^{+}$ controls molecule movement range and has a symmetric distribution. The movement size parameter ${r_{move}}$ affects the efficiency of sampling. If ${r_{move}}$ is too small, the acceptance ratio (which is the ratio of the number of proposed states from $q$ to the number of accepted states) becomes very high. However, in this case, the correlation between accepted states is strong. Thus independent states cannot be attained in an efficient manner. On the other hand, if ${r_{move}}$ is too large, independent states are easily sampled but proposal states are rarely accepted.

The procedure of standard MH is as follows:

Set sample steps $s=0$. 2. 2.

Set initial molecule positions ${\bm{r}}^{N}(s)$. 3. 3.

Calculate probability $p_{pre}=p({\bm{r}}^{N}(s=0))$. 4. 4.

Repeat the process process below until a sufficient number of states have been accepted. 5. 5.

Generate state ${\bm{r}}^{N}_{*}$ from the proposal distribution $q$ and compare to the previously accepted state ${\bm{r}}^{N}(s)$. 6. 6.

Calculate the probability $p_{*}=p({\bm{r}}^{N}_{*})$. 7. 7.

Accept ${\bm{r}}^{N}_{*}$ as ${\bm{r}}^{N}(s+1)$ with a probability of $\min\left(1,p_{*}/p_{pre}\right)$. 8. 8.

If accepted $s\leftarrow s+1$, $p_{pre}\leftarrow p_{*}$.

Through this process ${\bm{r}}^{N}(s)$ follows $p({\bm{r}}^{N})$.

Molecular Flow Metropolis-Hastings (MFMH)

The main idea behind the MFMH method is to avoid high energy states in the proposal distribution. We transform the proposal distribution $q({\bm{r}}_{i*}|{\bm{r}}_{i})$(Eq.14) into $q^{\prime}({\bm{r}}_{i*}|{\bm{r}}_{i})$ using the inverse square flow. Such proposal distribution $q^{\prime}$ generates no states which have extremely high energy, therefore, the acceptance ratio may be improved with smaller correlation. Furthermore, if the distribution $q$ is symmetric, the distribution of $q^{\prime}$ is also symmetric. This is because the instantaneous change is always space-constant through transformation and therefore the ratio of $q$ to $q^{\prime}$ is also constant. Thus, the calculation of the probability of $q^{\prime}$ is not required.

The procedure of MFMH is as follows:

Set sample steps $s=0$. 2. 2.

Set initial molecule positions before transform ${\bm{h}}^{N}(s)$. 3. 3.

Perform transformation into ${\bm{r}}^{N}(s=0)$ from ${\bm{h}}^{N}(s=0)$ using inverse square flow. 4. 4.

Calculate probability $p_{pre}=p({\bm{r}}^{N}(s=0))$. 5. 5.

Repeat the process process below until a sufficient number of states have been accepted. 6. 6.

Generate state ${\bm{h}}^{N}_{*}$ from the proposal distribution $q$ and compare to the previously accepted state ${\bm{h}}^{N}(s)$. 7. 7.

Perform transformation into ${\bm{r}}^{N}_{*}$ from ${\bm{h}}^{N}_{*}$ using inverse square flow. 8. 8.

Calculate the probability $p_{*}=p({\bm{r}}^{N}_{*})$. 9. 9.

Accept ${\bm{r}}^{N}_{*}$ as ${\bm{r}}^{N}(s+1)$ and ${\bm{h}}^{N}_{*}$ as ${\bm{h}}^{N}(s+1)$ with a probability of $\min\left(1,p_{*}/p_{pre}\right)$. 10. 10.

If accepted $s\leftarrow s+1$, $p_{pre}\leftarrow p_{*}$.

Note that, for MFMH, the integration time, $T$, for performing transformations must be same among all states. By following this process we can ensure that ${\bm{r}}^{N}(s)$ follows $p({\bm{r}}^{N})$.

5. Simulation results

Transformation of uniformly distributed molecules

In the first test of our method, we begin by verifying the results for transformations by the inverse square flow in the presence of periodic boundary conditions. As seen in Fig.1, systems that have several molecules pairs in close proximity are transformed into systems that exclude these molecule pairs. To demonstrate this clearly, we performed inverse square flow transformations on states of uniformly distributed molecules, and calculated the radial distribution functions (RDF), comparing the initial and transformed configurations. The RDF $g(r)$ describes the probability of finding a particle pair at a separation distance of $r$, relative to that of an ideal gas. The RDF $g(r)$ is defined mathematically as

[TABLE]

where $\left<n(r)\right>$ denotes the average number of particle-pairs with separation distances between $r$ and $r+dr$, $\rho$ is average density of the system. The configuration of molecules ${\bm{r}}^{N}$ in the system before the transformation was sampled from a uniform distribution. The results of these tests are shown in Fig.3. The results in Fig.3 indicates that the RDF profiles for the transformed states vary with integration time, $T$. The blue line ($T=0$) represents the RDF of initial configuration (uniformly distributed molecules). As time $T$ becomes longer, the separation distance between molecules becomes larger and molecule pairs in close proximity, which are the cause of extremely high energy states, are completely excluded. Therefore, these transformed states are suitable for efficient sampling.

Sampling Lennard-Jones particles system by MFMH method

Next, we confirmed the efficiency of the MFMH method using by simulation of a Lennard-Jones(LJ) particle system. This systems contains only monoatomic molecules whose interactions depend only on the LJ potential (?). For LJ particles, the two-body intermolecular potential energy $U_{ij}(r)$ is represented in dimensionless units as:

[TABLE]

where $r=\|{\bm{r}}_{i}-{\bm{r}}_{j}\|$. To test the efficiency of the MFMH method, we compared the acceptance ratios between the standard MH and the MFMH methods for the same proposal distribution $q$ (Eq.14) and molecule movement range ${r_{move}}$. The target system contains $N=64$ LJ particles, and the system $Temperature$ and particle density $\rho$ are 4 and 1, respectively. The results of these simulations are shown in Fig.4 and Fig.5. We generated $10^{6}$ proposal states using the standard MH method, and $10^{6}$ proposal states using the MFMH method. The molecule movement range was set to ${r_{move}}=0.05$. As a result, the energy distribution of states for both methods showed good agreement, except for the case where $T=0.160$. Further to this, the acceptance ratios were found to be 0.112 for the standard MH method and 0.244 for the MFMH method at $T=0.080$. Consistency of the energy distributions indicates that the MFMH method accurately sampled the system ensemble. These energy distributions were also consistent with previous research (?). The reason that the case for $T=0.160$ produces in an incorrect distribution can be considered a result of the fact that MF avoids excessively low energy states, which are necessary for correct MH sampling. Furthermore, even though the same value of ${r_{move}}$ was used, the acceptance ratio for the MFMH method is, in the best case, more than twice that of the standard MH. It may be claimed with some justification, therefore, that the MFMH method samples states more efficiently, and with the same level of accuracy, as the standard MH method.

6. Conclusions

In this work, we have proposed a novel and efficient sampling method “Molecular Flow Monte Carlo” using inverse square flow. Inverse square flow transforms any distribution into a new distribution which precludes the possibility of molecule pairs occurring in extremely close proximity, with a constant change in probability density. Therefore, we can avoid high energy states from the proposal distribution of the Metropolis-Hasting method, which improves efficiency. This improvement of efficiency has been confirmed through case study simulations of molecular systems.

Appendix

A. Non-derivative form of the pairwise force of inverse square flow

[TABLE]

may be rewritten as

[TABLE]

by performing differentiation.

B. Instantaneous change of the periodic sum of inverse square flow

We now investigate the instantaneous change of the following flow,

[TABLE]

${\bm{f}}^{p}({\bm{r}})$ is constructed from a real-space term and reciprocal-space term. The real-space term is

[TABLE]

Incidentally, if $f_{ij}(r){\bm{r}}$ is represented as ${\bm{\nabla}}_{\bm{r}}h_{ij}(r)$ by using some function $h_{ij}$, Eq.1 can be expressed as

[TABLE]

for $D=3$. Using Eq.22, and assuming ${\bm{f}}^{p}({\bm{r}})$ has only real-space terms, the instantaneous change may be calculated as

[TABLE]

by setting

[TABLE]

If we now think about the screening factor $G$ in the limit of $G\rightarrow 0$, the reciplocal-space term goes to [math], therefore

[TABLE]

holds. This equation holds any value of $G\in(0,\infty)$ because the screening factor $G$ only divides the inverse square flow into real and reciprocal terms and their sum is constant.