ここでは2次元の古典XY模型を考えることにします。
2次元XY模型は通常の相転移はありせんが、渦（トポロジー）と関係した Berezinskii-Kosterlitz-Thouless 相転移が存在します。
ちなみに通常の相転移が無いことには、スピンが連続値であることに伴なう Mermin-Wagnerの 定理が関係しています。
XY模型は量子色力学にも関係しますし、面白い応用としては音楽理論への適用もあります。
スピン模型は物理分野だけにとどまらず、経済、社会、生物など、その応用先は多岐に渡る非常に面白い模型です。
いかではc言語の表記を利用しながらその記述を見ていきましょう。
2次元 XY 模型のエネルギーは
E= -J \sum_<a,b> S_a S_b = -J \sum_<a,b> cos(θ_a - θ_b)
で与えられます。J は結合定数、Sは連続値を取るスピンで、θは対応する位相（0～2π）です。
ここで、\sum_<i,j> は最近接のスピンでのみ和を取ることを意味します。
具体的に座標を二つの整数 i,j で与えると、
E= -J \sum_{i,j} [ cos(θ_{i,j} - θ_{i-1,j}) + cos(θ_{i,j} - θ_{i+1,j}) + cos(θ_{i,j} - θ_{i,j-1}) cos(θ_{i,j} - θ_{i,j+1}) ]
と表記されます。
この模型の物理量の期待値を計算するにはどうすればよいでしょうか？
一番シンプルな方法は（計算時間はかかりますが）、イジング模型のページで説明しているメトロポリス法を用いる方法です。
イジング模型ではスピンが2値（+1と-1）でしたが今回は連続値ですので、「スピンを反転させる」を「0～2πの間での乱数を生成して足す」に変更する必要があります。
つまり、
spin = fmod(spin + (double)rand()%(2*PI), 2*PI)
の様に範囲を制限して 0～2πの範囲のランダムな数を足せばよいわけです。逆に言えば、それだけで計算できるようになります。
但しθは連続値を取るため、メトロポリス法でのスピンの更新回数はイジング模型に比べて格段に多く取らなけらばなりません。
2次元XY模型はBerezinskii-Kosterlitz-Thouless 相転移が存在するため、例えば比熱を計算しただけでは相転移温度を見積もることはできません。
その場合は相関関数 G(r-r')=<cos(θ_r - θ_r')> を計算し、そのデータを ∝exp(-R/ξ) の関数でフィットすることで相関長ξを数値的に見出します。
その後、相関長が発散する温度を計算すると相転移温度が（十分な注意の下に）分かります。
ちなみに通常の相転移では相関長が冪的に発散しますが 、Berezinskii-Kosterlitz-Thouless相転移では指数的に発散します。
そのため、exp(a/(T-b)) のような関数でフィットすると良いでしょう。ただし、有限体積効果には注意してください。
相関長ではなく helicity modulus という量からも相転移温度を調べることができます。これは捻りに対する応答をみる量となっています。
また、persisitent homology というトポロジカルデータ解析の手法を用いても相転移温度を見出すことができるという報告もあります。
色々な量を計算して、どの量がBerezinskii-Kosterlitz-Thouless 相転移に関係するのかを探るのもよい卒研（場合によっては修論）のテーマになると思います。
かなり適当に作ったファイルも公開しておきます。
＊ c言語で書いていますが仕様の違いのため Visual Studio 等ではそのままではコンパイルできないと思います（gccではコンパイル確認済み）。
＊非常に重いデータファイルが生成されるので、仕様を理解したうえで使ってください。
＊パラメータファイルは自分でいじって下さい。コードを読むとどうすればよいかわかります。
＊比熱はあまり精度良く計算できていません。参考論文に書かれているようにエネルギーを数値微分すると良いでしょう。