重力場中にあるブラウン粒子の平衡モンテカルロシミュレーション

某餃子

2024年9月27日 19:18

卒業研究に取り掛かる過程で練習した平衡モンテカルロ法によるシミュレーションについてまとめていきます。
内容としては初歩的です。

設定は以下の通りです。

2次元直交座標系を考え、$${y}$$軸(縦軸)下方向に一様な重力が作用する
系内には$${N}$$個のブラウン粒子がある
粒子間の相互作用は無視
系の温度$${T}$$は一定

平衡モンテカルロ法

Hoover氏の著書である計算統計力学を参照して、今回の設定を基にした平衡モンテカルロ法の概要を以下に述べます。

ランダムに粒子１つを選ぶ
選んだ粒子を適当に動かす
動かした後の粒子のポテンシャルエネルギー$${U_{new}}$$を計算し、動かす前のポテンシャルエネルギー$${U_{old}}$$との差分$${\Delta U = U_{old} - U_{new}}$$を取る
ボルツマン因子$${e^{-\beta\Delta U}}$$と、ある確率分布に由来する$${[0, 1)}$$で正規化された実数値を取る確率変数$${p}$$との大小関係を比較する。
4.において、もし$${e^{-\beta\Delta U}>p}$$ならば粒子の移動を受け入れ、$${e^{-\beta\Delta U}\leq p}$$ならば粒子は動かさない
上記の1.～5.までを$${N}$$個の全粒子に対して行い、十分な回数繰り返す

ただし、コンピュータでは無次元化した式を用いて上記の手順を踏むので、実際にプログラム中で使用する式は少し異なります。

次に、式の無次元化を行っていきます。

式の無次元化

計算の中で現れる定数と式を整理します。

粒子の質量：$${m}$$
重力加速度：$${g}$$
系の温度：$${T}$$
ボルツマン定数：$${k_{B}}$$
粒子のポテンシャルエネルギー：$${U=mgy}$$

このうち変数は$${y, U}$$ですので、これらを無次元化します。
無次元化にはそれぞれの次元にあった基準量が必要なので、作っていきます。
まずエネルギーの基準量ですが、ここは素直に$${k_{B}T}$$を採用します。
次に長さの次元を持つ物理量ですが、これも上記の定数の積で作っていきます。それぞれの次元を観察すると、$${\frac{k_{B}T}{mg}}$$が妥当そうです。よって、$${l_{ref}=\frac{k_{B}T}{mg}}$$として、これを長さの基準量とします。

これらを使って無次元化した変数を定義していきます。

$$
\begin{align*}
\tilde{U} &= \frac{U}{k_{B}T} = \frac{mgy}{k_{B}T} \\
\tilde{y} &= \frac{y}{l_{ref}} = \frac{mgy}{k_{B}T}
\end{align*}
$$

よって、無次元化後は$${\tilde{U} = \tilde{y}}$$となります。
ボルツマン因子内の$${\Delta U}$$と$${k_{B}T}$$の比の部分も、無次元化後は$${\Delta \tilde{U}}$$と$${1}$$の比となります。

上記の無次元化した式を使って、シミュレーションプログラムをpythonで実装します。

プログラム

計算に当たって、境界条件として周期境界条件を使いました。

import numpy as np
import matplotlib.pyplot as plt
import random

def emc_fbc(num_particles, particles, width, height, move_width):
    new_particles = particles.copy()
    keys = [i for i in range(num_particles)]
    random.shuffle(keys)

    for i in keys:
        dx = random.uniform(-1, 1) * move_width
        dy = random.uniform(-1, 1) * move_width

        new_x = particles[i, 0] + dx
        new_y = particles[i, 1] + dy


        # 境界条件(固定端)
        if new_x < 0:
            new_x += dx*2
        if new_x > width:
            new_x -= dx*2
        if new_y < 0:
            new_y += dy*2
        if new_y > height:
            new_y -= dy*2

        # エネルギー変化
        u_old = particles[i, 1]
        u_new = new_y
        delta_u = u_new - u_old

        # 受け入れ判断
        if  np.exp(-delta_u) > random.uniform(0, 1):
            new_particles[i, 0] = new_x
            new_particles[i, 1] = new_y

    return new_particles


# 粒子の数と空間の大きさ
num_particles = 1000
width = 10
height = 10

# 初期状態の設定（ランダムに配置）1列目がx座標, 2列目がy座標
x = np.random.rand(num_particles) * width
rng = np.random.default_rng()
y = rng.normal(height/2, 0.4, num_particles)
#y = np.random.rand(num_particles) * height
particles = np.column_stack((x, y))
init_x, init_y = particles[:, 0], particles[:, 1]

# シミュレーションのパラメータ
num_steps = 1000
move_width = 1.0

# シミュレーションの実行
for step in range(num_steps):
    particles = emc_fbc(num_particles, particles, width, height, move_width)

print('finish')

計算結果

上記のプログラムを走らせた結果を載せます。

次に縦軸を適当な長さで区切って、それぞれの区切りの中にある粒子の数を数えてみます。区切る長さは適当に$${ 10/9=1.111… }$$としました。

米沢富美子氏の著書「ブラウン運動」の72ページには、今回のシミュレーションと同様の系において「(中略)… ペランは同様の実験を数多く行って, いずれの場合にも濃度の比が等比級数の形になることを確かめた. 」とあります。今回の結果は程々にそれを再現できていると思います。

記事更新(2024/11/7) :
プログラム訂正(境界条件を周期境界から固定端に)，←に伴って結果のグラフ(1枚目)を変更，gif削除，分布をヒストグラムに変更