ワングとランダウのアルゴリズム

王富高とデイビッド・P・ランドウ^[1]が提案したワン・ランダウアルゴリズムは、システムの状態密度を推定するために設計されたモンテカルロ法です。この方法は、非マルコフランダムウォークを実行し、利用可能なすべてのエネルギースペクトルをすばやく訪問して状態密度を構築します。ワン・ランダウアルゴリズムは、マルチカノニカルシミュレーションを実行するために必要な状態密度を取得するための重要な方法です。

ワン・ランダウアルゴリズムは、コスト（またはエネルギー）関数によって特徴付けられるあらゆるシステムに適用できます。たとえば、数値積分の解法^[2]やタンパク質の折り畳みに適用されています。^{[3] [4]ワン・ランダウサンプリングは}メタダイナミクスアルゴリズム に関連しています。^[5]

ワングとランダウのアルゴリズムは、コスト関数によって特徴付けられるシステムの状態密度の推定値を取得するために使用されます。これは、漸近的にマルチカノニカルアンサンブルに収束する非マルコフ確率過程を使用します。^[1] (つまり、状態密度の逆のサンプリング分布を持つメトロポリス-ヘイスティングスアルゴリズム) 主な結果は、このサンプリング分布により、エネルギー障壁が見えないシミュレーションがもたらされることです。これは、アルゴリズムがメトロポリスアルゴリズムよりもはるかに速く、アクセス可能なすべての状態 (好ましい状態と好ましくない状態) を訪れることを意味します。^[6]

位相空間 上に定義されたシステムと、推定される状態密度 を持つスペクトル で制限されるコスト関数 E (たとえば エネルギー) について考えます。推定量は です。Wang と Landau のアルゴリズムは離散スペクトルで機能するため、^[1]スペクトルはの差を持つ N 個の離散値に分割され、${\displaystyle \オメガ}$${\displaystyle E\in \Gamma =[E_{\min },E_{\max }]}$${\displaystyle \rho (E)}$${\displaystyle {\hat {\rho }}(E)\equiv \exp(S(E))}$${\displaystyle \Gamma}$${\displaystyle \Delta }$

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

この離散スペクトルが与えられた場合、アルゴリズムは次のように初期化されます。

• ミクロカノニカルエントロピーのすべてのエントリをゼロに設定すると、${\displaystyle S(E_{i})=0\ \ i=1,2,...,N}$
• 初期化と${\displaystyle f=1}$
• ランダムな構成を設定することで、システムをランダムに初期化します。${\displaystyle {\boldsymbol {r}}\in \Omega }$

次に、アルゴリズムはマルチカノニカルアンサンブルシミュレーションを実行します。^[1]システムの位相空間におけるメトロポリス-ヘイスティングスランダムウォークで、確率分布は で与えられ、新しい状態を提案する確率は確率分布 で与えられます。訪問したエネルギーのヒストグラムが保存されます。メトロポリス-ヘイスティングスアルゴリズムと同様に、提案-承認ステップが実行され、次のようになります（メトロポリス-ヘイスティングスアルゴリズムの概要を参照）。 ${\displaystyle P({\boldsymbol {r}})=1/{\hat {\rho }}(E({\boldsymbol {r}}))=\exp(-S(E({\boldsymbol {r}})))}$${\displaystyle g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')}$${\displaystyle H(E)}$

1. 任意の提案分布に従って状態を提案する${\displaystyle {\boldsymbol {r}}'\in \Omega }$${\displaystyle g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')}$
2. 提案された状態を承認/拒否する

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

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

各提案承認ステップの後、システムはある値に遷移し、1 ずつ増加して次の更新が実行されます。 ${\displaystyle E_{i}}$${\displaystyle H(E_{i})}$

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

これはアルゴリズムの重要なステップであり、Wang と Landau のアルゴリズムを非マルコフ的にするものです。つまり、確率過程は過程の履歴に依存するようになります。したがって、次にその特定のエネルギーを持つ状態への提案があった場合、その提案は拒否される可能性が高くなります。この意味で、アルゴリズムはシステムにスペクトル全体を均等に訪問するように強制します。^[1]その結果、ヒストグラムはますます平坦になります。ただし、この平坦さは、計算されたエントロピーが正確なエントロピーにどれだけよく近似されているかに依存し、これは当然 f の値に依存します。^[7]正確なエントロピー (およびヒストグラムの平坦さ) をよりよく近似するために、M 回の提案-承認ステップの後に f が減らされます。 ${\displaystyle E_{i}}$${\displaystyle H(E)}$

${\displaystyle f\leftarrow f/2}$。

その後、fを常に2で割って更新すると飽和誤差が生じる可能性があることが示されました。^[7]この問題を回避するためにWangとLandauの方法に小さな修正を加えると、 に比例するf係数が使用されます。ここで、 はシミュレーションのステップ数に比例します。^[7]${\displaystyle 1/t}$${\displaystyle t}$

調和振動子ポテンシャルの DOS を取得したいと考えています。

${\displaystyle E(x)=x^{2},\,}$

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

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

最後の積分を実行すると、

${\displaystyle \rho (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の累乗で与えられ、指数はシステムの次元の関数になります。

したがって、状態密度の解析形式はすでにわかっているので、単純な調和振動子ポテンシャルを使用して Wang-Landau アルゴリズムの精度をテストできます。したがって、 Wang-Landau アルゴリズムによって得られた推定状態密度を と比較します。 ${\displaystyle {\hat {\rho }}(E)}$${\displaystyle \rho (E)}$

以下は、対称提案分布 g が使用されていることを前提とした、 Pythonでの Wang-Landau アルゴリズムのサンプルコードです。

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

コードでは、研究対象の基礎となるシステムである「システム」を考慮します。

currentEnergy  =  system.randomConfiguration () #ランダムな初期設定  

while  f  >  epsilon : 
    system.proposeConfiguration () # 提案された構成が提案されるproposalEnergy = system.proposedEnergy () #提案された構成のエネルギーが計算される  
        

    if  random ()  <  exp ( entropy [ currentEnergy ]  -  entropy [ proposalEnergy ]): 
        # 受け入れられた場合、エネルギーとシステムを更新します: 
        currentEnergy  =  proposalEnergy 
        system . acceptProposedConfiguration () 
    else : 
        # 拒否された場合
        system . rejectProposedConfiguration ()

    H [現在のエネルギー]  +=  1
    エントロピー[現在のエネルギー]  +=  f

    if  isFlat ( H ):   # isFlat はヒストグラムが平坦かどうかをテストします (例: 95% 平坦) 
        H [:]  =  0 
        f  *=  0.5   # f パラメータを調整します

分子動力学（MD）は通常モンテカルロ（MC）よりも好まれるため、フラットエネルギーサンプリングのための基本的なWLの考え方を組み込んだMDアルゴリズムが望ましい。そのアルゴリズムは、ボストン大学のJaegil Kimらによって 開発された統計温度分子動力学（STMD）である^{[8] 。}

重要な最初のステップは、統計的温度モンテカルロ (STMC) アルゴリズムによって達成されました。WLMC では、状態密度を直接操作するため、システム サイズに応じてエネルギー ビンの数が大幅に増加します。STMC は、統計的温度 という強度量を中心にしています。ここで、E は位置エネルギーです。関係 と組み合わせると(ここで と設定)、状態密度を更新するための WL ルールは、離散化された統計的温度を更新するためのルールになります。 ${\displaystyle T(E)=1/(dS(E)/dE)}$${\displaystyle \オメガ(E)=e^{S(E)}}$${\displaystyle k_{B}=1}$

${\displaystyle {\tilde {T}}_{j\pm 1}^{\prime }=\alpha _{j\pm 1}{\tilde {T}}_{j\pm 1},}$

ここで、はエネルギービンのサイズ、は実行中の推定値を表します。 ^[1]では、fを次のように定義します。システムがそのビンのエネルギーにアクセスしたときにi番目のエネルギービンのDOSの推定値を乗算する係数>1。 ${\displaystyle \alpha _{j\pm 1}=1/(1\mp \delta f{\tilde {T}}_{j\pm 1}),\delta f=(lnf/2\Delta E),\Delta E}$${\displaystyle {\tilde {T}}}$

詳細は文献^[8]に記載されています。 の初期推定値と範囲が から の間に制限されている場合、シミュレーションは WLMC の場合と同様に進行しますが、数値的には大きな違いがあります。 の補間により、その逆関数を積分して推定された の連続表現が得られ、WL よりも大きなエネルギー ビンを使用できます。 の異なる値は、受け入れ確率を評価するときに同じエネルギー ビン内で使用できます。 ヒストグラムの変動が平均の 20% 未満の場合、 は に従って減少します。 ${\displaystyle T(E)}$${\displaystyle T_{L}}$${\displaystyle T_{U}}$${\displaystyle {\tilde {T}}(E)}$${\displaystyle S(E)}$${\displaystyle S(E)}$${\displaystyle f}$${\displaystyle f\rightarrow {\sqrt {f}}}$

STMC は、イジング モデルとレナード ジョーンズ液体について WL と比較されました。エネルギー ビンのサイズが大きくなると、STMC はかなりの範囲で同じ結果を得ますが、WL のパフォーマンスは急速に低下します。STMD は、より迅速な収束のために、より小さな初期値を使用できます。つまり、STMC では、同じ品質の結果を得るために必要な手順が少なくなります。 ${\displaystyle f_{d}=f-1}$

ここで、主な結果である STMD について考えてみましょう。これは、がすべての位置を表す、ポテンシャル エネルギー から導出された力による温度での標準 MD シミュレーションにおいて、構成のサンプリング ウェイトは であるという観察に基づいています。さらに、力が関数 から導出されている場合、サンプリング ウェイトは です。 ${\displaystyle T_{0}}$${\displaystyle E([x])}$${\displaystyle [x]}$${\displaystyle e^{-E([x])/T_{0}}}$${\displaystyle W(E)}$${\displaystyle e^{-W(E([x]))/T_{0}}}$

フラットエネルギーサンプリングの場合、有効ポテンシャルを- エントロピー分子動力学とします。すると重みは です。状態密度は なので、それらの積はフラットエネルギーサンプリングを与えます。 ${\displaystyle T_{0}S(E)}$${\displaystyle e^{-S(E)}}$${\displaystyle e^{+S(E)}}$

力は次のように計算される。

${\displaystyle F=(-d/dx)T_{0}S(E)=T_{0}(dS/dE)(-d/dx)E([x])=(T_{0}/T(E))F^{0}}$

ここで、 はポテンシャル エネルギーから得られる通常の力を表します。 通常の力を係数でスケーリングすると、フラットなエネルギー サンプリングが生成されます。 ${\displaystyle F^{0}}$${\displaystyle (T_{0}/T(E))}$

STMDは、定数およびVで通常のMDアルゴリズムから開始します。力は示されているようにスケーリングされ、統計温度はSTMCと同じ手順を使用して、時間ステップごとに更新されます。シミュレーションがフラットエネルギーサンプリングに収束するにつれて、実行中の推定値は真のに収束します。収束を早める手順を含む技術的な詳細は、^[8]および^{[9]に記載されています。}${\displaystyle T_{0}}$${\displaystyle {\tilde {T}}(E)}$${\displaystyle T(E)}$

STMD では、通常どおり速度を制御するため運動温度と呼ばれますが、構成サンプリングには入りません。これは通常とは異なります。したがって、STMD は高速粒子で低エネルギーを調べることができます。任意の標準平均は再重み付けによって計算できますが、統計温度 は追加の分析なしですぐに利用できます。これは相転移の研究に非常に役立ちます。有限ナノシステムでは、あらゆる「サブフェーズ転移」に対応する特徴があります。十分に強い転移の場合、 の S ループ上の等面積構築によって転移温度が得られます。 ${\displaystyle T_{0}}$${\displaystyle T(E)}$${\displaystyle T(E)}$${\displaystyle 1/T(E)}$

STMDはBUグループによって改良され、^[9]彼らと他の研究者によっていくつかのシステムに適用されました。D. Stelterは、私たちが強度量の作業に重点を置いているにもかかわらず、は強度が強いことを認識しました。しかし、は強度が強く、ヒストグラムの平坦性に基づく手順は、固定数の時間ステップごとに半分にカットすることに置き換えられました。この単純な変更により、STMDは完全に強度が強くなり、大規模システムのパフォーマンスが大幅に向上します。^[9]さらに、強度の最終値は、収束した の誤差の大きさを決定する定数であり、システムサイズに依存しません。STMDは、LAMMPSではfix stmdとして実装されています。 ${\displaystyle ln(f)}$${\displaystyle \delta f=(ln(f)/2\Delta E)}$${\displaystyle f\rightarrow {\sqrt {f}}}$${\displaystyle \delta f}$${\displaystyle \delta f}$${\displaystyle T(E)}$

STMD は特に相転移に有効です。転移を起こすには過冷却または過熱が必要なので、標準的なシミュレーションでは平衡情報を得ることはできません。しかし、STMD を実行すると、低エネルギー状態または高エネルギー状態に閉じ込められることなく、加熱と冷却の自然な進行を伴うフラットなエネルギー サンプリングが得られます。最近では、脂質で包まれたナノ粒子の流体/ゲル転移^[9]に適用されています。

レプリカ交換STMD ^[10]もBUグループによって発表されている。