論文「Determining the density of states for classical statistical models」の実装
またまた友人からモンテカルロ法系の論文を教えてもらったので読んで実装。
その論文は
これは”Wang-Landauアルゴリズム”と呼ばれている計算手法で、統計物理学でいう”状態密度”、つまり”同じエネルギーを持つ状態の数はどのくらいあるのか？”を表す量を推計するために使用するアルゴリズム。正確には”密度”なんで数じゃないけど。状態密度を一旦計算してしまうと、いちいちモンテカルロを走らせ直さなくても異なる温度での物理量が簡単に計算しなおせて便利ですよとそういうご利益がある手法。
ちょいと調べるとこれはwikipediaにも独立な項目があるくらい有名なもんらしい。
さらにこれを多次元の数値積分を計算するのに応用している例も。
Analysis of the convergence of the 1/t and Wang-Landau algorithms in the calculation of multidimensional integrals
ここではまたオレオレ適当モデルを作成して実際に正しく計算できるかを試しに実装。
モデルは以下のように作成
- 全状態数は６
- 状態１〜３まではエネルギー３をとる（高エネルギー値）
- 状態４〜５まではエネルギー２をとる（中エネルギー値）
- 状態６はエネルギー１をとる（低エネルギー値）
- 従って
- エネルギー３をとる状態は３つ（状態密度単位×３）
- エネルギー２をとる状態は２つ（状態密度単位×２）
- エネルギー１をとる状態は１つ（状態密度単位・・・と設定）
となる。
状態密度を推定した結果はこちら
左から順に「低エネルギーの状態密度」・「中エネルギーの状態密度」・「高エネルギーの状態密度」*1となっており、実際に設定した状態密度（比）とほぼ一致することがわかる。
以下、コード。
# -*- coding: utf-8 -*- import numpy import matplotlib.pyplot #対象とするモデルの定義 class System: def __init__(self): #現在と一個前の系の状態 self.state_current = 1 self.state_previous = 1 #系のエネルギーが何個（パターン）あるか self.number_of_energies = 3 #乱数生成機 self.rand = numpy.random.mtrand.RandomState(100) #系の状態に合わせてエネルギー用インデックスとそのエネルギー値自身を返却 def GetEnergy(self, state): if state in [1, 2, 3]: return (2, 3) elif state in [4, 5]: return (1, 2) else : return (0, 1) def GetCurrentEnergy(self): return self.GetEnergy(self.state_current) #系の状態をランダムに変化させる def GoNextState(self): self.state_previous = self.state_current self.state_current = self.rand.random_integers(1,6) #系の状態を一個前のに戻す def BackPreviousState(self): self.state_current = self.state_previous #ヒストグラムが平滑化の判定 def IsFlat(histgram): average = numpy.average(histgram) if all(map(lambda z : abs(z - average) / average < 0.05, histgram)): return True else: return False def EstimateDOS(system, eps_log_modification_factor): #乱数ジェネレーター rand = numpy.random.mtrand.RandomState(100) #各エネルギーに対してインデックスを持つヒストグラムと状態密度 histgram = [0 for i in range(system.number_of_energies)] log_density_of_state = [0 for i in range(system.number_of_energies)] #LOG版の修正ファクター？とその誤差の許容範囲の設定 log_modification_factor = 1 #モンテカルロステップカウンター。定期的にヒストグラムの平滑化具合をチェックするのに使用 counter_mc = 0 while log_modification_factor > eps_log_modification_factor: counter_mc += 1 #状態変化前後でのエネルギー値を取得 index_energy_current, _ = system.GetCurrentEnergy() system.GoNextState() index_energy_new , _ = system.GetCurrentEnergy() log_dos_new = log_density_of_state[index_energy_new] log_dos_current = log_density_of_state[index_energy_current] if rand.rand() < min(1.0, numpy.exp(log_dos_current - log_dos_new)): log_density_of_state[index_energy_new] += log_modification_factor histgram[index_energy_new] += 1 else: system.BackPreviousState() log_density_of_state[index_energy_current] += log_modification_factor histgram[index_energy_current] += 1 if(counter_mc % 1000 == 0): counter_mc = 0 if IsFlat(histgram): histgram = [0 for i in range(system.number_of_energies)] log_modification_factor *= 0.5 return log_density_of_state if __name__ == '__main__': #システム生成 x = System() #収束計算の許容範囲 eps_log_modification_factor = 10**(-8) #状態密度計算 log_density_of_state = EstimateDOS(x, eps_log_modification_factor) density_of_state = [i / numpy.exp(log_density_of_state)[0] for i in numpy.exp(log_density_of_state)] print density_of_state #plot index = [i + 1 for i in range(x.number_of_energies)] width = 0.35 pos = matplotlib.pyplot.bar(index, density_of_state, width, color = 'r', align = 'center') matplotlib.pyplot.ylabel('Ratio(normalized)') matplotlib.pyplot.title('Density of states') matplotlib.pyplot.show()
*1:低エネルギーの状態密度で規格化
