（様式 20）成果報告書

TSUBAME 共同利用 平成 24 年度 学術利用 成果報告書

利用課題名 TSUBAME2 GPU によるスピン系のクラスターアルゴリズム・モンテカルロシミュレーション
英文： Cluster algorithm Monte Carlo simulations of spin systems using TSUBAME2 GPU
利用課題責任者 岡部 豊
First name Surname Yutaka Okabe

Affiliation

所属 首都大学東京理工学研究科
Graduate School of Science and Engineering, Tokyo Metropolitan University

邦文抄録（300 字程度）
格子スピンモデルは、物質内の相転移現象を解析するのに使われる。格子スピンモデルの解析にはよくマルコフ連鎖モンテカルロ法が用いられており、今日までに多くのマルコフ連鎖モンテカルロ法が提案されている。近年、GPU による計算の高速化が多くの分野で注目されている。本プロジェクトではマルコフ連鎖モンテカルロ法の一つ、Swendsen-Wang (SW) マルチクラスターアルゴリズムの複数の GPU を用いたアルゴリズムの開発を行う。本年度は昨年度作成したプログラムを拡張し、格子スピンモデルの一種、2 次元 XY モデルに対して、大規模かつ高速計算を行った。結果、従来問題とされていた有限サイズ補正の問題を解決し、転移温度に対しても従来よりも一桁精度高い値の算出に成功した。
英文抄録（100 words 程度）
The classical spin model is one model for understanding mechanisms of phase transition. The classicalspin model is often analyzed by using a Markov Chain Monte Carlo simulation, and many algorithms
for Markov Chain Monte Carlo simulation have been proposed. In this paper we study thetwo-dimensional classical XY model by the large-scale Monte Carlo simulation of the Swendsen-Wang
multi-cluster algorithm using multiple GPUs on the open science supercomputer TSUBAME 2.0.
Simulating systems up to the linear system size L = 65536, we investigate the Kosterlitz–Thouless(KT) transition. The obtained inverse KT temperature KT is 1.11996(6). We estimate the exponent tospecify the multiplicative logarithmic correction, -2r, and precisely reproduce the theoretical prediction-2r =1/8

Keywords: Monte Carlo Simulation, classical spin model

背景と目的

は昨年度作成したプログラムを拡張させ、実際に問題多体系の物理的性質を理解するためにモンテカルロとなっている系での計算に取り組むことを目指した。具法が標準的な手法として広く用いられているが、シング体的に、格子スピンモデルの一種、2 次元 XY モデルにルスピンフリップのモンテカルロ法には転移温度付近で対し、大規模かつ高速なシミュレーションを行うことで今緩和時間が急激に大きくなる問題があり、これを解消まで難問とされていた問題を解決した[2]。

する一つの手法としてクラスターアルゴリズムが提唱さ2 次元 XY モデルは Kosterlitz–Thouless (KT)転移れている。しかし、クラスターアルゴリズムの大規模な[3,4]という独特な相転移を示す。通常の相転移では転計算には様々な問題があり、今までの報告ではクラス移温度上で相関長が発散するが、KT 転移では転移温ターアルゴリズムは大規模化な計算に適していないと度以下の全ての温度で相関長が発散する。通常の相されていた。本プロジェクトでは、格子スピンモデルに転移と比較して、KT 転移での転移温度付近での相関対し、複数の GPU を使った大規模な Swedesn-Wang

長は速く発散し(SW) クラスターアルゴリズム・モンテカルロシミュレーションを開発し、クラスターアルゴリズムの大規模化の c  c   A exp   t t (T )  exp (1)

基礎の確立と計算の高速化の２つを目的とする。昨年
度はアルゴリズムの考案及び実装を行った[1]。本年度のように振る舞う。ここで t=(T-TKT)/ TKT、c は定数、

（様式 20）成果報告書
TKT は KT 転移温度である。一方、2 次元 XY モデルのることが導かれており、数値計算では以下のように定モンテカルロシミュレーション解析には有限サイズによ義される[10,11]。

る対数補正から生じる困難も予想されている[4,5]。KT e 

転移温度上での磁化率は

 m2L2 L L2(ln L)

2 r(2)と見積もられる。ここでの L は一辺の格子サイズであり、JL2 2sT(4)e1 cos( i   j )L2 i , j (5)e1 cos( i   j )L2 i , j (6)

η=1/4、r=-1/16 と理論的に予測されている[4]。モンテ

ここでの式(4)での < > は温度平均を意味する。このカルロシミュレーションを用いた数値解析も多く試みら基準から各サイズの KT 転移温度を求め、有限サイズれているが、対数補正項 r の算出は未だに成功していスケーリングより、補正項 r の値および KT 転移温度にない。対数補正項 r の歴史は参考文献[6]に示されていついての計算を行う。計算を行った格子サイズは L=64,る。モンテカルロシミュレーションの解析が成功していな

128, 256, 512, 1024, 2048, 4096, 8192, 16384,

い理由としては計算サイズが小さいことが最大の理由

32768, 65536 であり、L=8192 以上の格子サイズの計である。今まで大きいサイズの計算でも L=2048 までの算に複数の GPU を使用した。各サイズの磁化率と KT計算しかない[7]。そのため、申請者らは 2 次元 XY モ

転移温度の値を図 1 に示す。

デルに対して、複数の GPU を用いた L=655536 までの大規模かつ高速計算を行い、この難問の解決を図る

L

βKT(L)<m>L103 MCSことを試みた。

64

1.09167(25)0.34924(29)

4000

128

1.09788(32)0.29695(36)

4000

256

1.10224(21)0.25245(28)

4000

512

1.10530(26)0.21427(27)

4000

1024

1.10768(19)0.18195(16)

4000

2048

1.10966(20)0.15443(26)

4000

4096

1.11097(08)0.13085(15)

4000

8192

1.11211(27)0.11086(32)

2000

16384

1.11314(20)0.09395(25)1000

め、補正項 r の値を議論するためには KT 転移温度上32768

1.11389(21)0.07955(22)500でのシミュレーションが必要となる。しかし、KT 転移温

65536

1.11457(14)0.06737(14)400

結果および考察二次元 XY モデルのハミルトニアンを以下に示す。

H   J  cos( i   j )(3)

i , J 

ここでのθi はサイト i での XY スピンの角度を表し、0から 2πの値をとる。和は全ての最近接格子で行う。格子としては正方格子をとり、周期的境界条件をとる。
磁化率の式(2)は KT 転移温度上での振る舞いのた度の値については定まっておらず、昨年度開発した

図 1: 各サイズの KT 点転移温度と磁化率

SW クラスターアルゴリズムをそのまま適用出来ない。
そのため、本研究では SW クラスターアルゴリズムを拡

図 1 の括弧は誤差を表している。また、有限サイズスケ張させた確率変動クラスターアルゴリズム[8]を用いてーリングを用いた各サイズの KT 転移温度の関係式を計算を行う。確率変動クラスターアルゴリズムはある基式(7)に、また各サイズの補正項の関係式を式(8)にそ準を定め、その基準に対して温度を変化させるアルゴ

れぞれ示し、それぞれの様子を図 2,3 に示す。

リ ズ ム で あ る 。 今 回 、 そ の 基 準 と し て helicity

c 2  KT KT ( L)   KT (ln bL) 2modulus[4,9]を用いる。helicity modulus は繰り込み群方程式から、KT 転移温度上で 2/πから 0 に変化す(7)

（様式 20）成果報告書

ln( m 2L

L1 / 4 )  const.  2r ln(ln aL)(8)

式(7)の b は定数であり、βは温度の逆数である。

従来では莫大な時間がかかった計算を高速に計算することが出来た。今後はこのような有限サイズ補正のあるモデルに対しての計算を行なっていく。また、本年度はこの研究に加え、パーコレーション問題への応用に1.12ついても取り組んだ。
1.11996−0.826(ln bL)−2

これらの結果については、TSUBAME ESJ Vol.7の英文および和文の解説記事で報告した。

1.11

KT(L)

なお、本研究は大学院生の小村幸浩君との共同研究である。
1.1

参考文献

1.09
00.01

0.02

−20.03(ln bL)[1] Komura, Y. and Okabe, Y. (2013). Comp. Phys.
Comm. 184, 40.

図 2: 各サイズの KT 転移温度の有限サイズスケーリ[2] Komura, Y. and Okabe, Y. (2012). J. Phys. Soc.
Jpn. 81, 113001.

ング[3] Kosterlitz, J. M. and Thouless, D. J. (1973). J.
Phys. C: Solid State Physics 6, 1181.

0.08

21/4ln (<m >LL )[4] Kosterlitz, J. M. (1974). J. Phys. C: Solid State
Physics 7, 1046.

0.06[5] Janke, W. (1997). Phys. Rev. B 55, 3580.
[6] Kenna, R. (2006).

0.04

Cond. Mat. Phys. 9, 283.

[7] Hasenbusch, M. (2005).
−0.262+0.128 ln(ln(aL))

0.02

J. Phys. A: Math.

Gen. 38, 5869.
[8] Tomita, Y. and Okabe, Y. (2002).

022.2

2.4

2.6ln(ln(aL))
図 3: 各サイズの磁化率の有限サイズスケーリング

Phys. Rev. B

66, 180401(R).
[9] Nelson D. R. and Kosterlitz, J. M. (1977). Phys.
Rev. Lett. 39, 1201.
[10] Teitel, S. and Jayaprakash, C. (1983).

図 2 より、KT 転移温度 1/ TKT が 1.11996(6)と求ま

Rev. B 27, 598.

り、従来の値、1.1199(1) [7]、1.1200(1) [12]よりも一[11] Weber, H. and Minnhagen, P. (1988).

桁精度高い値の算出に成功した。また、図 3 より対数

Rev. B 37, 5986.

補正項 r=-1/16 になることが確かめられた。

[12] Arisue, H. (2009).

まとめ、今後の課題
本年度は昨年度作成した複数の GPU を用いた大規模な SW クラスターアルゴリズムを拡張させ、実際に問題となっていた 2 次元 XY モデルの計算を行った。このモデルの問題を解決するためには大規模かつ高速な計算が不可欠であり、複数の GPU を使用することでPhys.
Phys.

Phys. Rev. E 79, 011107.