ワン・ランダウ法

ワン・ランダウ法（ワン・ランダウほう、英: Wang and Landau algorithm）は、Fugao Wangとデビッド・P・ランダウにより発案された、系の状態密度を計算するために用いられるモンテカルロ法のひとつである。このアルゴリズムでは状態密度の計算に必要な、系の取り得る全ての状態のエネルギーを迅速に計算するため、非マルコフ連鎖ランダムウォークを行う。マルチカノニカル法の実行に必要となる、状態密度を計算するために重要である。

ワンランダウ法はコスト（エネルギー）関数により特徴づけられるどのような系にも応用可能である。たとえば、数値積分やタンパク質フォールディングへの応用が知られている。ワン・ランダウサンプリングはメタダイナミクス法とも関連する。

概要

ワン・ランダウ法により、コスト関数で特徴づけられる系の状態密度を得ることができる。非マルコフ連鎖確率過程を用いて、漸近的にマルチカノニカルアンサンブルを得る（サンプル分布が状態密度の逆数となるようなメトロポリス・ヘイスティングス法に漸近する）。このことの帰結として、エネルギー障壁を無視したシミュレーションを行うことができる。つまり、このアルゴリズムは（コスト関数的に良いものも悪いものも含めた）系の取り得る全ての状態をメトロポリス法よりも非常に速く網羅することができるのである。

アルゴリズム

位相空間 $\Omega$ 上に定義される系と、そのコスト関数（例えばエネルギー） E を考える。 E のスペクトルは $E\in \Gamma =[E_{\min },E_{\max }]$ の範囲に限定されるものとし、付随する状態密度 $\rho (E)\equiv \exp(S(E))$ を計算する。 ワン・ランダウ法は離散スペクトルを扱うため、スペクトル $\Gamma$ をそれぞれ $\Delta$ だけ異なる N 個の離散値で分割する。

$N={\frac {E_{\max }-E_{\min }}{\Delta }}$

このような離散スペクトルに対し、まず次の初期化を行う。

• 全ての区間のミクロカノニカルエントロピーをゼロとする。

  $S(E_{i})=0\ \ i=1,2,...,N$

• f = 1 とする。
• 系の初期配置 ${\boldsymbol {r}}\in \Omega$ をランダムに選ぶ。

次に、マルチカノニカルアンサンブルシミュレーションを行う。つまり、位相空間上でランダムウォークを行い、確率分布

$P({\boldsymbol {r}})=1/\rho (E({\boldsymbol {r}}))=\exp(-S(E({\boldsymbol {r}})))$

を再現するような遷移確率 $g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')$ でメトロポリスシミューレーションを行う。このとき、訪れた状態のエネルギーをヒストグラム $H(E)$ に記録する。そしてメトロポリス法と同様に新たな配置の生成と採択/棄却プロセスを実行する。

1. 新たな配置 ${\boldsymbol {r}}'\in \Omega$ を任意の遷移確率分布 $g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')$ に従って生成する。
   
2. 新たな配置を次の確率で採択/棄却する。

   $A({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')=\min \left(1,e^{S-S'}{\frac {g({\boldsymbol {r}}'\rightarrow {\boldsymbol {r}})}{g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')}}\right)$

   ここで、 $S=S(E({\boldsymbol {r}}))$ および $S'=S(E({\boldsymbol {r}}'))$ とする。

このステップを終えたとき、系の遷移先が $E_{i}$ だったとすると $H(E_{i})$ を 1 だけ増やし、エントロピーを以下のように更新する。

$S(E_{i})\leftarrow S(E_{i})+f$

このステップが本アルゴリズムの中心である。また、このステップが存在するため確率過程が過程の履歴に依存し、ワン・ランダウ法は非マルコフ連鎖となる。したがって、このエネルギー $E_{i}$ を持つ状態が次に生成された際にはより棄却されやすくなる。このようにして、本アルゴリズムはスペクトル全域を等しく訪れるよう強制する。この帰結として、ヒストグラム $H(E)$ はシミュレーションを進めるにつれてより平坦となる。しかし、この平坦さは計算されたエントロピーがどれだけ実際のエントロピーを近似できているかに依存する。そして、近似の良さは f に依存する。実際のエントロピーに近づくため（すなわちヒストグラムを平坦にするため）に、 f を M ステップごとに次のように減らす。

$f\leftarrow f/2$

後に、 f を常に半分にし続けると飽和誤差を招くことが明らかとなった。この問題を避けるため、t をシミュレーションの総ステップ数として、 f を 1/t にする変更がアルゴリズムに加えられた。

アルゴリズムのテスト

今、次のような調和振動子ポテンシャルの DOS を求めたいとする。

$E(x)=x^{2}$

解析的には DOS は次のように与えられる。

$g(E)=\int \delta (E(x)-E)\,dx=\int \delta (x^{2}-E)\,dx$

最後の積分を解くと位相空間の次元によって次のような結果が得られる。

$g(E)\propto {\begin{cases}E^{-1/2},{\text{for }}x\in \mathbb {R} ^{1}\\{\text{const}},{\text{for }}x\in \mathbb {R} ^{2}\\E^{1/2},{\text{for }}x\in \mathbb {R} ^{3}\\\end{cases}}$

一般的に、多次元調和振動子の DOS は E のべき乗で与えられ、その指数は系の次元の関数となる。

このように単純な調和振動子ポテンシャルに対する状態密度は既に分かっているので、このポテンシャルに対してワン・ランダウ法を行い、 得られた状態密度 $\rho (E)$ と $g(E)$ を比べることによりワン・ランダウ法の精度を確かめることができる。

サンプルコード

次に示すのはPythonで実装されたワン・ランダウ法のサンプルコードである。遷移確率密度は対称であることを仮定している。

$g({\boldsymbol {x}}'\rightarrow {\boldsymbol {x}})=g({\boldsymbol {x}}\rightarrow {\boldsymbol {x}}')$

このコードにおいて、対象となる系は "system" により表されるものとしている。

currentEnergy = system.randomConfiguration() # a random initial configuration

while (f > epsilon):
    system.proposeConfiguration() # a proposed configuration is proposed
    proposedEnergy = system.proposedEnergy() # the energy of the proposed configuration computed

    if (random() < exp(entropy[currentEnergy]-entropy[proposedEnergy])):
        # if accepted, update the energy and the system:
        currentEnergy = proposedEnergy
        system.acceptProposedConfiguration()
    else:
        # if rejected
        system.rejectProposedConfiguration()
    
    H[currentEnergy] += 1
    entropy[currentEnergy] += f
    
    if (isFlat(H)): # isFlat tests whether the histogram is flat (e.g. 95% flatness)
        H[:] = 0
        f *= 0.5 # refine the f parameter

ワン・ランダウ分子動力学法

ワン・ランダウ法はモンテカルロ法のみならず分子動力学法にも適用することができる。そのためには、系の温度を次のように制御する。

$T'(E)\rightarrow \left({\frac {\partial S(E)}{\partial E}}\right)T(E)$

ここで S(E) は系のエントロピー、 T(E) はミクロカノニカル温度、 T′(E) は分子動力学シミュレーションに実際用いられる「スケールされた」温度である。