Difference between revisions of "Tutorials:Code-01 Python/ja"

From ALPS
Jump to: navigation, search
m
m (Text replace - "http://alps.comp-phys.org/static/tutorials2.1.0" to "http://alps.comp-phys.org/static/tutorials2.2.0")
 
Line 5: Line 5:
 
このチュートリアルでは、Python-ALPS を用いてどのようにしてシミュレーションコードを記述するかを学びます。計算物理における"hello world"として、古典2次元イジングモデルの局所更新モンテカルロ法をとりあげます。
 
このチュートリアルでは、Python-ALPS を用いてどのようにしてシミュレーションコードを記述するかを学びます。計算物理における"hello world"として、古典2次元イジングモデルの局所更新モンテカルロ法をとりあげます。
  
モンテカルロシミュレーションの典型的な構造を記述したスケルトンコード([http://alps.comp-phys.org/static/tutorials2.1.0/code-01-python/ising-skeleton.py ising-skeleton.py])を元にしながら、一歩一歩完全なコードを作っていきます。
+
モンテカルロシミュレーションの典型的な構造を記述したスケルトンコード([http://alps.comp-phys.org/static/tutorials2.2.0/code-01-python/ising-skeleton.py ising-skeleton.py])を元にしながら、一歩一歩完全なコードを作っていきます。
  
 
最初に、必要なPythonモジュールをインポートします。
 
最初に、必要なPythonモジュールをインポートします。

Latest revision as of 22:09, 28 September 2013

Pythonを利用したALPS入門

このチュートリアルでは、Python-ALPS を用いてどのようにしてシミュレーションコードを記述するかを学びます。計算物理における"hello world"として、古典2次元イジングモデルの局所更新モンテカルロ法をとりあげます。

モンテカルロシミュレーションの典型的な構造を記述したスケルトンコード(ising-skeleton.py)を元にしながら、一歩一歩完全なコードを作っていきます。

最初に、必要なPythonモジュールをインポートします。

import math
import pyalps
import pyalps.alea as alpsalea
import pyalps.pytools as alpstools

まず、Simulationクラスを実装します。これは、シミュレーションの初期化をするメソッド、シミュレーションを行うメソッド、結果をHDF5ファイルに格納するメソッドなどが含まれます。

class Simulation:
# Seed random number generator: self.rng() will give a random float from the interval [0,1)
   rng = alpstools.rng(42)
   
   def __init__(self,beta,L):
       self.L = L
       self.beta = beta
       
       # Init exponential map
       self.exp_table = dict()
       for E in range(-4,5,2): 
         self.exp_table[E] = math.exp(2*beta*E)
       
       # Init random spin configuration
       self.spins = [ [2*self.randint(2)-1 for j in range(L)] for i in range(L) ]
       
       # Init observables
       self.energy = alpsalea.RealObservable('E')
       self.magnetization = alpsalea.RealObservable('m')
       self.abs_magnetization = alpsalea.RealObservable('|m|')

__init__メソッドはSimulation<クラスのオブジェクトがどのようにインスタンス化されるのかを定義します。引数として格子サイズLと逆温度betaが渡されます。これらのパラメータに基づいて取りうるボルツマン重み(の比)を計算しておき、また、ランダムな配置で正方格子上のイジングスピンを初期化します。さらに求める物理量の初期化をおこないます。ここではALPS aleaライブラリが物理量の解析をおこなうようにpython-ALPSフレームワークを使用しています。クラス内では乱数生成器の初期化と乱数の種の設定もおこないます。乱数生成器にはメルセンヌ・ツイスタMT19937を使用しています。これは十分長い周期と良い統計的性質を持っており、モンテカルロ法に使うのに適しています。

def save(self, filename):
       pyalps.save_parameters(filename, {'L':self.L, 'BETA':self.beta, 'SWEEPS':self.n, 'THERMALIZATION':self.ntherm})
       self.abs_magnetization.save(filename)
       self.energy.save(filename)
       self.magnetization.save(filename)

saveメソッドはシミュレーションパラメータと結果をHDF5ファイルに書き出します。

def run(self,ntherm,n):
       # Thermalize for ntherm steps
       self.n = n
       self.ntherm = ntherm
       while ntherm > 0:
           self.step()
           ntherm = ntherm-1
       # Run n steps
       while n > 0:
           self.step()
           self.measure()
           n = n-1
       # Print observables
       print '|m|:\t', self.abs_magnetization.mean, '+-', self.abs_magnetization.error, ',\t tau =', self.abs_magnetization.tau
       print 'E:\t', self.energy.mean, '+-', self.energy.error, ',\t tau =', self.energy.tau
       print 'm:\t', self.magnetization.mean, '+-', self.magnetization.error, ',\t tau =', self.magnetization.tau

runメソッドはstepメソッドを呼び出して更新を行い、またmeasureメソッドを呼び出して物理量の測定を行います。熱平衡化中は測定を行わず、更新のみをします。すべての計算が完了したら、得られた物理量の平均値、統計誤差、自己相関時間を出力します。

def step(self):
    for s in range(self.L*self.L):
      # Pick random site k=(i,j)
       ...
      # Measure local energy e = -s_k * sum_{l nn k} s_l
       ...		
      # Flip s_k with probability exp(2 beta e)
      ...

モンテカルロ掃引はstepメソッドでおこなわれます。メトロポリスアルゴリズムでは、スピンをランダムに選択し、確率p_{accept} = min(1,e^{-\beta \Delta E})で反転します。\Delta Eは初期配位と提案配位でのエネルギー差です。この手順はL^2回繰り返されます。メトロポリスアルゴリムの実装は良い練習になるでしょう。次のように定義されたrandint関数が使えます。

def randint(self,max):
       return int(max*self.rng())

物理量の測定はmeasureメソッドに実装します。

def measure(self):
   E = 0.	# energy
   M = 0.	# magnetization
   for i in range(self.L):
       for j in range(self.L):
           E -= ...
           M += ...
   # Add sample to observables
   self.energy << E/(self.L*self.L)
   self.magnetization << M/(self.L*self.L)
   self.abs_magnetization << abs(M)/(self.L*self.L)

エネルギーと磁化をスピン配置から計算し、ALPS observableに加えます。実装してみてください。

メトロポリス更新や物理量の測定のコードを実装したら、vispythonalpspythonを使って計算を実行することができます。今回は、\beta = 1/k_B Tの異なる値でのスキャンをおこないます。

L = 4	# Linear lattice size
N = 5000	# of simulation steps
print '# L:', L, 'N:', N
# Scan beta range [0,1] in steps of 0.1
for beta in [0.,.1,.2,.3,.4,.5,.6,.7,.8,.9,1.]:
   print '-----------'
   print 'beta =', beta
   sim = Simulation(beta,L)
   sim.run(N/2,N)
   sim.save('ising.'+str(beta)+'.h5')


python-ALPSを使用する利点は、もし、U = <A>/<B>のような物理量を評価する場合、ジャックナイフ法が実行されて、自動的に正しい平均値や誤差を得られることにあります。たとえば、今回のコードを拡張し、m^2m^4を測定量として加え、これらからBinderキュムラントU_4=< m^4> /< m^2>^2を計算します。通常、Binder キュムラントは有限サイズスケーリングを用いて臨界温度を求めるのに使用するので、新たにL=6,8でも計算します。

これら2つの系での計算が正常に実行されたら、まず最初にデータをロードし、各システムサイズLごとに、逆温度\beta の関数として m^2m^4を得ます。

data = pyalps.loadMeasurements(pyalps.getResultFiles(pattern='ising.L*'),['m^2', 'm^4'])
m2=pyalps.collectXY(data,x='BETA',y='m^2',foreach=['L'])
m4=pyalps.collectXY(data,x='BETA',y='m^4',foreach=['L'])

Binder キュムラントの計算をし、datasetのリストに格納します。

u=[]
for i in range(len(m2)):
    d = pyalps.DataSet()
    d.propsylabel='U4'
    d.props = m2[i].props
    d.x= m2[i].x
    d.y = m4[i].y/m2[i].y/m2[i].y
    u.append(d)

次のコマンドで、Binder キュムラントをプロットできます。

import pyalps.plot as plt 

plt.figure()
plt.plot(u)
plt.xlabel('Inverse Temperature $\beta$')
plt.ylabel('Binder Cumulant U4 $g$')
plt.title('2D Ising model')
plt.show()


どんなグラフがあらわれましたか?もしさらに相転移や有限サイズスケーリングについて学びたいのであれば、 MC-07 Phase transition in the Ising model を参照してください。