2018-02-02

ASE: 窒素分子のCu(111)の二原子層における吸着構造の最適化

Python module のASE を使って、吸着構造を最適化。
ここでは簡単な有効媒質理論を使うが、calculatorを選ぶことで、AbinitとかVASPとかも使えるらしい。
Introduction: Nitrogen on copper — ASE documentation
元サイトのチュートリアルに少し追加して、吸着した時に結合距離がどう変わるかを見る。
吸着エネルギーを出しているが、ここではあまり関係ない。

from ase import Atoms
from ase.calculators.emt import EMT
from ase.constraints import FixAtoms
from ase.optimize import QuasiNewton
from ase.build import fcc111, add_adsorbate
from ase.visualize import view

h = 1.85 # initial adsorption distance
d = 1.10 # initial bond distance

# make Cu structure
slab = fcc111('Cu', size=(4, 4, 2), vacuum=10.0)
view( slab )

# calculate Cu energy
slab.set_calculator(EMT())
e_slab = slab.get_potential_energy()

# make N2 structure
molecule = Atoms('2N', positions=[(0., 0., 0.), (0., 0., d)])
molecule.set_calculator(EMT())
##e_N2 = molecule.get_potential_energy()

# optimise N2 structure
dyn = QuasiNewton(molecule, trajectory='N2.traj')
dyn.run(fmax=0.05)
print( molecule.get_all_distances() )

# calculate N2 energy
e_N2 = molecule.get_potential_energy()

# make a whole structure
add_adsorbate(slab, molecule, h, 'ontop')
view( slab )
print( [ dis for dis in slab.get_all_distances()[-2] if dis < 3. ] ) #  distance Cu-N1, N1-N1, N1-N2

# optimize adsorption
constraint = FixAtoms(mask=[a.symbol != 'N' for a in slab])
slab.set_constraint(constraint)
dyn = QuasiNewton(slab, trajectory='N2Cu.traj')
dyn.run(fmax=0.05)

print('Adsorption energy:', e_slab + e_N2 - slab.get_potential_energy())
view( slab )
print( [ dis for dis in slab.get_all_distances()[-2] if dis < 3. ] ) #  distance Cu-N1, N1-N1, N1-N2

add_adsorbateで足された原子群はリストの後ろから入る。
距離を出力する部分は、ちょっとサボって、3オングストローム以下だけ出力するようにした。
本当は、短い距離順にソートして、そこから三つ持ってくるのが良いと思う。

結果：
[1.8499999999999996, 0.0, 0.99809885325797332]
[2.0610952626456882, 0.0, 1.0480140019810147]

1.85は適当に決めた吸着距離だからどうでもよい。
N-Nの距離に注目すると、吸着前は 0.998 だったが、吸着後は 1.048 に伸びることが予言された。
「Cuに近い方の窒素原子の方がより引っ張られるから、結合が伸びる」というイメージを支持するものと思われる。

今はN2分子がCu表面に対して垂直だが、横にした時どうかとか、そういう比較をし始めるときには、Adsorption energyが必要になってくるだろう。

2018-01-13

Path operator と Coherent Potential Approximation (CPA) （誤り）

量子力学数学

#このページは誤っています。

Path operator を次のように用意しておく。

$\downarrow$ 誤り（一般的なPath operatorの定義ではない。）
$\displaystyle \tau = G + GtG + GtGtG + ... \\ \displaystyle = G \left( 1 - tG \right)^{-1} = \left( G^{-1} - t \right)^{-1} \\ \displaystyle = t^{-1} \left[ \left( 1 - tG \right)^{-1} - 1 \right]= \left( t^{-1} - G \right)^{-1} - t^{-1}$
$\uparrow$ 誤り

$G$ はサイト間の移動、 $t$ はサイト上での繰り込まれた相互作用（site T matrix）を表すとする。
後のために $t$ を一応、site potential $v$ を使って定義しておく。
$\displaystyle t = v + v g v + v + v g v g v + ... = \left( v^{-1} - g \right)^{-1}$
$g$ はサイト内の伝搬を表すと定義するが、後の議論には出てこないのでここでは割愛。

サイトと $t$ が一対一対応すれば良いが、例えば合金のように一つのサイトに対して複数の原子がある確率を持って占有する場合、 $t$ をsite potentialの平均から作る近似（仮想結晶近似:VCA）は粗過ぎる。
何故なら、本来は固定された $t$ のセットを全ての可能な配置に対して求め、それを配置数で割るという平均の取り方をするべきであり、複数種類の $t$ が同時にそのサイトに重みを持って存在している訳ではない。

したがって、配置に対して、すなわちpathの取り方に対して、ある種の平均操作を行うのが妥当だと思われる。
そのような平均操作によって得られたpath operatorを $\bar{ \tau }$ と定義する。
仮に $\bar{ \tau }$ が求めまっていると仮定する。
この $\bar{ \tau }$ を与えるような $\bar{ t }$ を上でpath operatorを用意したように、次の様に定義する。

$\downarrow$ 誤り（一般的なPath operatorの定義に沿っていない）
$\displaystyle \bar{ \tau } = G + G \bar{ t } G + ... = \left( G^{-1} - \bar{ t } \right)^{-1} \\ \displaystyle \therefore \; \bar{ t } = G^{-1} - \bar{ \tau }^{-1} \\ \displaystyle or \\ \displaystyle \bar{ t }^{-1} = \left( G^{-1} - \bar{ \tau }^{-1} \right)^{-1} = G G^{-1} \left( G^{-1} - \bar{ \tau }^{-1} \right)^{-1} \bar{ \tau }^{-1} \bar{ \tau } = G \left( \bar{ \tau } - G \right)^{-1} \bar{ \tau }$
$\uparrow$ 誤り

仮に求まっているとした平均的な場（と呼んで良いかは微妙だが、ここでは有効場と呼ぶことにする）に、注目している相互作用 $t_{\alpha}$ が埋まっているとする。
その場合のpath operator $\bar{ \tau }_{\alpha}$ は次の様に求められる*（厳密になんでそうなるかは不勉強）。
$\downarrow$ 誤り
$\displaystyle \bar{ \tau }_{\alpha} = \bar{\tau} \left[ 1 + \left( \bar{ t }^{-1} - t^{-1}_{\alpha} \right) \bar{\tau} \right]$
$\uparrow$ 誤り
$\downarrow$ 訂正
$\displaystyle \bar{ \tau }^{nn}_{\alpha} = \bar{\tau}^{nn} \left[ 1 + \left( t^{-1}_{\alpha} - \bar{ t }^{-1} \right) \bar{\tau}^{nn} \right] ^{-1} = \left[ \left( \bar{\tau}^{nn} \right)^{-1} + \left( t^{-1}_{\alpha} - \bar{ t }^{-1} \right) \right] ^{-1}$
$n$ はサイト番号を表す。 $\tau \neq \tau^{nn}$ に注意（前者は行列、後者はその要素）。
$\uparrow$ 訂正

これの意味は、次の変形を踏まえると分かり易くなる。
$\downarrow$ 誤り（一般的なPath operatorの定義に沿っていない）
$\displaystyle \bar{ t }^{-1} - t^{-1}_{\alpha} = \left( \bar{ v }^{-1} - g \right) - \left( v^{-1}_{\alpha} - g \right) = \bar{ v }^{-1} - v^{-1}_{\alpha} = \frac{ v_{ \alpha } - \bar{ v } }{ v_{ \alpha } \bar{ v } } = \frac{ \Delta v }{ v_{ \alpha } \bar{ v } }$
$\uparrow$ 誤り

つまり、以下の様に、有効場を与える様なsite potential（決してただの平均ではない（と思う））を考慮したいsite potentialで置き換えることで、「埋まる」ことを表現している（なんで $v_{\alpha} \bar{ v }$ で割るのかは不明。しかもこれだとただの一次摂動に思える）。
$\downarrow$ 誤り（一般的なPath operatorの定義に沿っていない）
$\displaystyle \bar{ \tau }_{\alpha} = \bar{\tau} + \bar{\tau} \frac{ \Delta v }{ v_{ \alpha } \bar{ v } } \bar{\tau}$
$\uparrow$ 誤り

有効場に埋まっている各相互作用、ここでは二元系として $A, B$ 、に対して、有効path operator $\tau_A, \tau_B$ が求まったが、これらが有効場を構成しているわけだから、この平均がもともとの有効場に一致すると考えれば、
$\displaystyle x_A \bar{ \tau }_A + x_B \bar{ \tau }_B = \bar{ \tau }$
$x_A, x_B$ は $A, B$ の割合である。

これで変数の閉じた自己無撞着方程式が得られた。
この方法でpath operatorを得る方法をCoherent Potential Approximation (CPA)と呼ぶ。

~~やっぱり、 $\bar{ \tau }_{\alpha}$ の導出を理解しないと、わかった気になれんなぁ。。。~~ Path operatorの定義が間違っていたからわからんのは当たり前。

*Ref: V. Popescu et al., J. Synchrotron Rad. 6 (1999) 711.

2018-01-11

グランドカノニカル平均からのKeldysh形式に移る前の導入の不満

量子力学

参考文献："Nonequilibrium Green's Functions Approach to Inhomogeneous Systems", Karsten Blazer and Michael Bonitz (Springer 2013).
該当箇所：2.1.1 Keldysh Contour

この手のもので自分が最初に詰まってしまうのは、

ハミルトニアンが時間に依存しているのかどうか

時間に依存する場合、時間の異なるハミルトニアンの交換を仮定しているのかどうか

が曖昧な点である。

最終的には、一般化された順序演算子を導入して好き勝手に順番を弄れるようにするから何でも良いのだろうが、導入の部分だけを真面目に考えると、適用範囲が狭過ぎて「これ意味あるの？？」と思ってしまった。でも導入だからあまり真面目に書かれていない気がする。

時間 $t$ におけるグランドカノニカル平均を次の様に定義。
$\displaystyle < \hat{A} > (t) = \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ e^{ - \beta ( \hat{H}(t) - \mu \hat{N} ) } \hat{A}_S \right\} = \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ e^{ - \beta \hat{K}( t ) } \hat{A}_S \right\}, \\ \displaystyle Z_0 = {\rm Tr } \left\{ e^{ - \beta \hat{K}( t_0 ) } \right\} \\ \displaystyle \hat{K}( t ) = \hat{H}( t ) - \mu \hat{N}$
ハットは演算子であることを表す（ハットの位置がズレて見えるのが不満で仕方がない）。

一般的な時間発展演算子は時間順序演算子 $\hat{T}$ を使って表すと、
${ \displaystyle \hat{U}( t_1, t_0 ) = \hat{T} {\rm exp} \left( - \frac{ i }{ \hbar } \int^{t_1}_{t_0} dt \; \hat{H}( t ) \right) }$

ハミルトニアンが時間に依存するが、 $\left[ \hat{H}(t_1), \hat{H}(t_2) \right] = 0$ の時（時間毎にエネルギーがwell defined）、
${ \displaystyle \hat{U}( t_1, t_0 ) = {\rm exp} \left( - \frac{ i }{ \hbar } \int^{t_1}_{t_0} dt \; \hat{H}( t ) \right) }$

最後に、ハミルトニアンが時間に依存しない場合、
${ \displaystyle \hat{U}( t_1, t_0 ) = {\rm exp} \left( - \frac{ i ( t_1 - t_0 ) }{ \hbar } \hat{H} \right) }$

ここで、 $t_0$ における温度因子を次の様に定義。
${\displaystyle \hat{U}( t_0 - i \hbar \beta, t_0 ) = {\rm exp} \left( - \beta \hat{K}( t_0 ) \right) }$
ここの取り扱いは $\hat{K}( t ) = \hat{K}( t_0 )$ として、ハミルトニアン（もしくは $\hat{K}$ ）が時間に依存しないとした時の時間発展演算子に対応させたものになる。

$\hat{U}$ の中のハミルトニアンを $\hat{K}$ に読み替え、ハミルトニアン（ $\hat{K}$ ）が時間に依存しなければ、
${\displaystyle \hat{U}( t_0 - i \hbar \beta, t_0 ) = \hat{U}( t, t_0 ) \hat{U}( t_0, t ) \hat{U}( t_0 - i \hbar \beta, t_0 ) = \hat{U}( t, t_0 ) \hat{U}( t_0 - i \hbar \beta, t_0 ) \hat{U}( t_0, t ) }$

これによって、
$\displaystyle < \hat{A} > (t) = \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ \hat{U}( t, t_0 ) \hat{U}( t_0 - i \hbar \beta, t_0 ) \hat{U}( t_0, t ) \hat{A}_S \right\} \\ \displaystyle = \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ \hat{U}( t_0 - i \hbar \beta, t_0 ) \hat{U}( t_0, t ) \hat{A}_S \hat{U}( t, t_0 ) \right\}$

これは $t_0 \rightarrow t \rightarrow t_0 \rightarrow t_0 - i \hbar \beta$ という時間発展を考えましょうということを表している。
でもこれ自身は、ハミルトニアンが時間に依存していない時の話だから、こんな仰々しいことしないで松原先生にお願いすることになると思われる。

多分大事なのは、もっと一般のハミルトニアンの時には、直接こうには書けないけれども、「後で順番整理するからとりあえず好きなように並べていいよ！」として同じような時間経路を組んだらどう計算出来るか、という流れなのだと思う。
まあ時間に依存しなきゃ順番なんて何でも良いから、その時に順序演算子は要らなくて、厳密にその形でかけるようになるのは当たり前と言えば当たり前ですが。。。

Formulationは触ったことありますが、実際に運用したことはないので、参考文献を読み進める次第です。

2018-01-10

調和振動子と平面波の比較

量子力学

零点振動している調和振動子と同じエネルギーを持つ平面波とで波動関数を比較する。
$\displaystyle l_0 = \sqrt{ \frac{ \hbar }{ m \omega } } \\ \displaystyle \xi = \frac{ x }{ l_0 } \\ \displaystyle \psi_n( x ) = A_n H_n\left( \frac{ x }{ l_0 } \right) e^{ - \frac{1}{2} \left( \frac{ x }{ l_0 } \right)^2 } = A_n H_n( \xi ) e^{ - \frac{ \xi^2 }{ 2 } } \\ \displaystyle A_n = \sqrt{ \frac{ 1 }{ 2^n \; n! } } \; \sqrt{ \frac{ l_0 }{ \sqrt{ \pi } } } \\ \displaystyle H_0 = 1 \\ \displaystyle E_0 = \frac{ \hbar \omega }{ 2 }$
最後ら辺、サボりました。

ここで、 $E_0 = E_k$ として自由平面波のエネルギーを代入して、波数と有効距離 $l_0$ を結びつけると、
$\displaystyle l_0 = \sqrt{ \frac{ \hbar }{ m } \frac{ \hbar }{ 2 E_k } } = \sqrt{ \frac{ \hbar^2 }{ m } \frac{ 2 m }{ 2 \hbar^2 k^2_0 } } = \frac{ 1 }{ k_0 }$
綺麗に有効距離の逆数になる。

これで求めた $k_0$ を持つ平面波（ $cos( k_0 x )$ ）を $n=0$ の調和振動子波動関数と比較する（規格化定数は無視）と、全体で一周期分をちょっと超えたかなくらいの波が調和振動子の中に入っていることがわかる。
f:id:koideforest:20180110005333p:plain

もう少し厳密に見てみると、 $k_0$ を逆に固定して、 $k_0 x_0 = \pi$ となる $x_0$ における調和振動子(係数は無視)の値を求めると、
$\displaystyle e^{ - \frac{1}{2} \left( \frac{ x }{ l_0 } \right)^2 } = e^{ - \frac{1}{2} ( k_0 x )^2 } \rightarrow e^{ - \frac{1}{2} ( k_0 x_0 )^2 } = e^{ - \frac{\pi^2}{2} } = 0.00719... \sim 0.01 \ll 1$
したがって、十分小さいと見なせるところまでの範囲で、平面波一周期分が入っていることがわかる。

思っていたよりも調和振動子と平面波の対応関係がしっかりしていてビックリした。
原点近傍でポテンシャルがゼロに近くなるから、その近傍で自由解に近づくから当然といえば当然であるが。
自由解だけ波動関数があるような領域は、運動エネルギー（波動関数）ではなくポテンシャル項によってエネルギーが補填されているんだなぁとしみじみ感じた。

2017-12-14

PythonのLaguerre陪関数について

プログラミング数学量子力学

Wikipediaには日頃からよくお世話になっているが、自分で確認するのは大事だなと感じた事件。

水素原子の波動関数の可視化とかカッコイイなぁと思い、ネットサーフィン中に以下のサイト様を訪問。
scipyの特殊関数 - 篠突く雨の日記

自分も早速コピーしてやってみて、ふむふむと散布図の便利さに舌鼓を打っていた。
ちょっとアレンジして、波動関数の正負で色を変えるようにしたら気付いてしまった。
「ん？1s波動関数に節があんぞ？？」
波動関数の部分はscipy.special.assoc_laguerreで実装されていて、原因はこいつであることがわかった。
何か引数がおかしいのかと思い、Wikipediaで水素原子の解析解におけるLaguerre陪関数を確認。
水素原子におけるシュレーディンガー方程式の解 - Wikipedia
$\displaystyle \rho \frac{ d^{ 2 } u( \rho ) }{ d\rho^{ 2 } } + ( 2l + 2 - \rho ) \frac{ d u( \rho ) }{ dx } + ( n - l - 1 ) u( \rho ) = 0 \\ \displaystyle \rho \frac{ d^{ 2 } u( \rho ) }{ d\rho^{ 2 } } + ( ( 2l + 1 ) + 1 - \rho ) \frac{ d u( \rho ) }{ dx } + ( ( n + l ) - ( 2l + 1 ) ) u( \rho ) = 0 \\ \displaystyle u( \rho ) = c L^{ 2 l + 1 }_{ n + l }( \rho ) \\ \displaystyle L^m_k( \rho ) = \frac{ d^m }{ d\rho^m } e^{ \rho } \frac{ d^k }{ d\rho^k } ( e^{ -\rho } \rho^k ) \\$
Legendreではlとmは直接打ち込めば良かったけど、Laguerreはそういえばややこしかったんだよなぁ、と遠い記憶が蘇ってきた。
しかし、可視化スクリプトの引数はこれと対応したものが入っている。
実は上付きと下付きは逆の表記になっているとかあるか？？と思い、他も参照。

ラゲールの陪多項式 - Wikipedia
$\displaystyle \left( x \frac{ d^{ k + 2 } }{ dx^{ k + 2 } } + ( k + 1 - x ) \frac{ d^{ k + 1 } }{ dx^{ k + 1 } } + ( n - k ) \frac{ d^k }{ dx^k } \right) L^k_n( x ) = 0 \\ \displaystyle L^k_n( x ) = \frac{ d^{ k } }{ dx^{ k } } L_n( x ) = \frac{ d^{ k } }{ dx^{ k } } \left( e^x \frac{ d^{ n } }{ dx^{ n } } ( x^n e^{ -x } ) \right)$
上がkだったり下がkだったりしてめちゃくちゃややこしいが、とりあえずLaguerre多項式（下付き添字のみ）を上付き添字の分だけ微分しまくるのが陪関数であることがわかった。
その意味で、水素原子におけるシュレーディンガー方程式の解 - Wikipediaの上付き下付きの意味はラゲールの陪多項式 - Wikipediaと同じ。

scipy.special.assoc_laguerreのサイトに行ってみた。
scipy.special.assoc_laguerre — SciPy v1.1.0 Reference Guide
見る限り、スクリプトの引数の順番も間違っていなさそうだった。

しかし、

scipy.special.assoc_laguerre( x, n + l, 2 * l + 1 )

はn=1, l=0で $-x + 1$ のプロットを示し、意味不明。

そして英語版のwikipediaに辿り着く。
Laguerre polynomials - Wikipedia
$\displaystyle xy'' + ( \alpha + 1 - x ) y' + \beta y = 0 \\ \displaystyle L^{ ( \alpha ) }_0( x ) = 1 \\ \displaystyle L^{ ( \alpha ) }_1( x ) = 1 + \alpha - x$
これや、、、これやったんや。。。

つまり、scipy.special.assoc_laguerreを水素様原子の波動関数に使うには、
$\displaystyle L^{ \alpha }_{ \beta } = L^{ 2 l + 1 }_{ n - l - 1 }$
ということになる。
確かにこれで1s (n=1, l=0)の時 $L^1_0 = 1$ となり、めでたく波動関数の節が無くなる。

お手軽に特殊関数が使える様になったとはいえ、油断は禁物ですな。

2017-12-14

モンテカルロ積分で水素原子のエネルギー期待値を確認

プログラミング量子力学

「スレーター型軌道の数値積分は大変」と呪文のように聞かされて続けて心にハードルが出来てしまったので、ここらで打破することを試みた。
前回、多次元のモンテカルロ積分をまとめた。
koideforest.hatenadiary.com
今回はスクリプトを組んで、どれくらいちゃんと動くのか見てみることにした。
目標は、非相対論水素原子の解析解を使って、ハミルトニアンの期待値を出す。
以下、スクリプトのポイントをまとめる。

運動エネルギーを出すところの微分は、中心差分法を使った。簡単便利。

積分範囲は、水素原子波動関数の値が0.01を下回る距離をrmaxと定義し、その範囲で水素原子の8分の1が入る箱を考えた。対称性から、その箱の中の積分値を8倍すれば全体の積分値が求まる。

そんなことをしたのは、乱数の範囲が $[ 0, 1 ]$ の0始まりなので、全体を平行移動したりして水素原子全体を箱の中に無理矢理入れるのは余り得策ではなく感じたから。今回は素直に水素原子の原点を0にした。

箱の中の全体で積分値を求めるから、乱数点に関する制限は必要ない。下手にif分を入れると遅くなるらしいので、そこまで一般化してスクリプトを作るのは余り良くないかもしれない。今回はそこの一般化は行っていない。

以下ソース。

import numpy as np
import time

def MonteCalro3d_ws( func, scales ): # ws: whole space
    M = 10**6
    vec_rand = np.random.rand( M * 3 ).reshape( M, 3 )
    # Caution: the number of samplings is still M, not M * 3
    vec_rand = np.array( scales ) * vec_rand
    tot = np.sum( [ func( v[0], v[1], v[2] ) for v in vec_rand ] )
    return np.prod( scales ) * tot / float( M )

def h1s( x, y, z ):
    vec = np.array( [ x, y, z] )
    r   = np.linalg.norm( vec )
    return 2. * np.exp( - r ) / np.sqrt( 4. * np.pi )

def d2func_dx2( func, x ):
    h = 1.0e-7
    return ( func( x + 2.*h ) - 2.*func(x) + func( x - 2.*h ) ) / (2.*h)**2

def integrand( x, y, z ):
    v = np.array( [ x, y, z] )
    r = np.linalg.norm( v )
    deriv2x = d2func_dx2( lambda x: h1s( x,     v[1], v[2] ), v[0] )
    deriv2y = d2func_dx2( lambda y: h1s( v[0],     y, v[2] ), v[1] )
    deriv2z = d2func_dx2( lambda z: h1s( v[0],  v[1],    z ), v[2] )
    ke = - ( deriv2x + deriv2y + deriv2z ) / 2.
    pot_n = - h1s( v2[0], v2[1], v2[2] ) / r
    return h1s( v[0], v[1], v[2] ) * ( ke + pot_n )

rmax = max( [ x_ for x_ in np.linspace( 0, 5, 1000 ) \
              if h1s( x_, 0., 0. ) / h1s( 0., 0., 0. ) > 0.01 ] )
scales = [ rmax, rmax, rmax ]

start = time.time()
result = 8. * MonteCalro3d_ws( integrand, scales )
print( "Hydrogen atom energy = {}".format( result )  )
end = time.time()
print( end - start )

結果：
-0.500215806919
148.051234961

原子単位系は原子単位系でも、Hartree unitで計算しているので、理論値は-0.5でちゃんと計算できている。

ちなみにHartreeかRydbergかは運動エネルギーが $1/2$ されているかどうかで判断できる。

$1/2$ の無いRydbergの場合には、核ポテンシャルの項も二倍して理論値は-1になる。

細かい値を忘れたが、ざっくり 1 Hartree ~ 27 eV, 1 Rydberg ~ 13 eV, 1 Hartree = 2 Rydberg である。最後の等式は厳密。

値の安定性はサンプリング数 $M$ に依存するが、 $M = 10^6$ であれば有効数字一桁は安心できると思う。

cythonとか使えばもっと速くなる。

2017-12-13

N次元のモンテカルロ積分

数学

一般的なアイデアとして、積分を体積 x 密度（割合）に焼き直すところがスタート。
$\displaystyle I = \int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) \ = \left( \int_{ \Omega } d \textbf{ r } \right) \ \left( \frac{ \int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) }{ \int_{ \Omega } d \textbf{ r } } \right) \ = V_{\Omega} < f >_{\Omega}$

密度（割合）は単位量当たりにどれだけ対象となる事物があるかを表し、いわゆる期待値として焼き直せる。連続量の期待値なので、（確率）密度 or 重みをかけた積分に対応。それを離散量で近似しましょうというノリ。このあたりは本当は測度論とかが顔を出すんだろうなぁと妄想していますが、不勉強なので回避。
$\displaystyle < f >_{\Omega} = \int_{ \Omega } d \textbf{ r } \; \rho_{ \Omega }( \textbf{ r } ) f( \textbf{ r } ) \ \sim \sum^{ m }_{ \mu_1 } \cdots \sum^{ m }_{ \mu_N } p_{\Omega}( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (N) }_{ \mu_N } ) f( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (N) }_{ \mu_N } )$

$m$ はサンプリング数またはメッシュの細かさ。トータルのメッシュ数は $M \equiv m^N$ で定義する。
一次元と三次元に関して先ほどの和を解して書いてみると、
$\displaystyle \sum^{ m }_{ \mu } p_{\Omega}( x_{ \mu } ) f( x_{ \mu } ) = p_{\Omega}( x_1 ) f( x_1 ) + p_{\Omega}( x_2 ) f( x_2 ) + \cdots + p_{\Omega}( x_m ) f( x_m ) \\ \displaystyle \sum^{ m }_{ \mu_x } \sum^{ m }_{ \mu_y } \sum^{ m }_{ \mu_z } p_{\Omega}( x_{ \mu_x }, y_{ \mu_y }, z_{ \mu_z } ) f( x_{ \mu_x }, y_{ \mu_y }, z_{ \mu_z } ) \\ \displaystyle \qquad = p_{\Omega}( x_1, y_1, z_1 ) f( x_1, y_1, z_1 ) + p_{\Omega}( x_2, y_1, z_1 ) f( x_2, y_1, z_1 ) \\ \displaystyle \qquad + \cdots + p_{\Omega}( x_m, y_m, z_{ m-1 } ) f( x_m, y_m, z_{ m-1 } ) + p_{\Omega} ( x_m, y_m, z_m ) f( x_m, y_m, z_m )$

一次元において、 $0 \leq x \leq a$ における積分を、等間隔のメッシュを切った長方形で近似する場合で書いてみれば、もう少し期待値の描像が明確になる。
$\displaystyle \int^{a}_{0} dx \; f(x) \sim \sum^{ m }_{ i = 1 } \frac{ a }{ m } f( ( i - 1 ) a / m ) \\ \displaystyle \int^{a}_{0} dx \; f(x) = a < f >_a \\ \displaystyle \therefore \; < f >_a \sim \sum^{ m }_{ i = 1 } \frac{ 1 }{ m } f( ( i - 1 ) a / m ) \ = \sum^{ m }_{ i = 1 } p f( a ( i - 1 ) / m ) \\ \displaystyle \qquad \sim \sum^{ m }_{ i = 1 } p f( a R_i )$

全空間において一様に離散化・サンプリングすると、それは平均値を取りましょうということになる。
そして、これの近似として、 $0 < R_i < 1$ で規格化された昇順の乱数列 $\{ R \}$ を定義して置き換えたものが、積分領域に制約のない場合のモンテカルロ積分となる。
近似関係が分かりやすいようにあえて昇順と言ったが、結局は和なので順番は関係ない。
とりあえず多重積分したいという場合にはこっちの考え方になるだろう。いちいちメッシュの切り方を自分で考えなくて良いので、手軽で良い。
$p$ が単純にサンプリング数の逆数なので、モンテカルロ積分の例題で有名な $\pi$ を求める問題（円の中に入った点の数を数える）は、この考えでは出来ない。

これまでの話は領域 $\Omega$ の中で閉じていた。
なので、今までの $p_{\Omega}$ はあまり確率の意味はなく、単純に重みの役割が大きかった。
この領域 $\Omega$ を扱い易い任意の箱 $\Box$ に拡張するが、その際に余分な領域での値がゼロになるような関数 $\Theta$ を定義する。
$\displaystyle \int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) \ = \int_{ \Box } d \textbf{ r } \; \Theta( \textbf{ r } ) f( \textbf{ r } ) = V_{ \Box } < \Theta f >_{ \Box } \\ \Theta( \textbf{ r } ) = \left\{ \ \begin{array}{ ll } 1 & ( \textbf{ r } \in \Omega ) \\ 0 & (otherwise) \end{array} \right.$

これを先ほどと同様に密度関数で表したのちに離散化すると、
$\displaystyle < \Theta f >_{\Box} = \int_{ \Box } d \textbf{ r } \; \rho_{ \Box }( \textbf{ r } ) \Theta( \textbf{ r } ) f( \textbf{ r } ) \ = \int_{\Box} d \textbf{ r } \; \rho( \textbf{ r } ) f( \textbf{ r } ) \\ \displaystyle \sim \sum^{ m }_{ \mu_1 } \cdots \sum^{ m }_{ \mu_N } p( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (1) }_{ \mu_N } ) f( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (N) }_{ \mu_N } )$

上と同様に一次元の例で考えれば、
$\displaystyle \int^{a}_{0} dx \; f(x) = \int^{b}_{0} dx \; \theta( a - x ) f(x) \sim b \sum^{ m }_{ i = 1 } \frac{ \sigma( ( i - 1 ) b / m ) }{ m } f( ( i - 1 ) b / m ) \\ \displaystyle \int^{a}_{0} dx \; f(x) = b < \theta_a f >_b \\ \displaystyle \therefore \; < \theta_a f >_b \sim \sum^{ m }_{ i = 1 } \frac{ \sigma( ( i - 1 ) b / m ) }{ m } f( ( i - 1 ) b / m ) \ = \sum^{ m }_{ i = 1 } p( b ( i - 1 ) / m ) f( b ( i - 1 ) / m ) \\ \displaystyle \qquad \sim \sum^{ m }_{ i = 1 } p( b R_i ) \; f( b R_i ) \\$

ここで、
$\displaystyle \sigma( x ) = \left\{ \begin{array}{ll} 1 & ( x \in [ 0, \; a ] ) \\ 0 & (otherwise) \end{array} \right.$
と定義した。

$\sigma$ 及び $p$ で、積分領域の調整を乱数点が領域に入った入らないかで行うことが出来る。積分領域が複雑な時には、この方法が威力を発揮することになる。
乱数自体は $[ 0, 1 ]$ の間でしか値が出ないでの、計算したい系に合わせてスケール因子を掛ける必要がある。
一つの次元を $c$ 倍する度に結果を $c$ 倍すれば正しい積分結果が得られる。

nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。