● 2次元イジング模型

「計算物理」の教科書でFORTRANで設計されている2次元イジング模型を C++言語で設計してみよう。イジング模型と熱浴法など理論については 旧版を参照のこと。ここでは新しい内容のみを紹介する。


■ プログラムの実装

プログラムを設計し始めよう。ソースファイルは ising.cc である。まずは諸定数の定義である。

#include "cip.h"
#include "nxgraph.h"

//---- definition of fundermental constants
//---- physical constants
const double J      = 1.0;           // bonding constant
const double H      = 0.0;           // external magnetic field
const int SIZE      = 16;            // size of lattice

//---- experimental settings
// settings about temperature
const double T_INIT = 0.50;          // initial temperature
const double T_FINA = 5.00;          // final   temperature
const double T_STEP = 0.25;          // increments

// settings about sweep the lattice
const int RELAX     = 400;           // steps till relaxation
const int INTERVAL  = 10;            // steps per one measure

// settings about number of data to measure fluctuations of ...
const int SAMPLES   = 40;            // energy, magnetization
const int GROUPS    = 20;            // susceptibility, specific heat

// file names of output data
const char* fnames[4] = {  "energy.dat", "magnet.dat", "suscep.dat", "spheat.dat" };
// titles of each quantities
const char* titles[4] = {  "Energy", "Magnetization", "Susceptibility", "Specific heat" };
// file pointers
static FILE* fptrs[4];

//---- index indicators
const int ENERGY    = 0;
const int MAGNET    = 1;
const int SUSCEP    = 2;
const int SPHEAT    = 3;

//---- graphics settings
const int    WIN_WIDTH  = 256;
const int    WIN_HEIGHT = 256;

SIZE が格子の一辺のスピンの数である。 プログラムでは各温度に設定後 RELAX 回だけ格子を操引して 平衡状態にする。その後は INTERVAL 回の操引でデータを採取する。 エネルギーと磁化が簡単に求まるが、実験は一回の測定で終る ものではない。SAMPLES 回の測定を重ねてエネルギーと磁化の 平均値と分散を求める。この分散から帯磁率と比熱が求まる。 この4つの量の平均値とエラーを得るため、さらにこの測定作業を GROUPS 回行う。

このプログラムでは統計量のデータ収集が重要である。 データの平均値(mean)、分散(variance)、標準偏差(standard deviation)、 そして平均値の偏差(error)を機能的に計算する機構が必要である。 そこで次のclass Statistic を作成する。

//---- declaration and definition of class Statistic
class Statistic
{
  int n;                // current total of data
  double sum1, sum2;    // current sum of data and its square
public:
  inline Statistic( void ){
    n = 0;    sum1 = sum2 = 0.0;
  }
  //-- adds a new datum
  inline void Add( double datum ){
    n++;    sum1 += datum;    sum2 += datum*datum;
  }
  //-- current mean value
  inline double Mean( void ){
    return( sum1/n );
  }
  //-- current variance value
  inline double Variance( void ){
    return( (sum2 - sqr(sum1)/n)/(n-1) );
  }
  //-- current deviation value
  inline double Deviation( void ){
    return( sqrt( Variance() ) );
  }
  //-- current deviation of mean value, namely error value
  inline double Error( void ){
    return( sqrt( Variance()/n ) );
  }
};

内容は単純である。Constructorに総和の初期化を行わせる。 Add() methodにデータをひとつづつ与えて、統計量が 必要な所でそのmethodを呼び出すだけである。この各種統計量の 定義式を以下にまとめておく。

次にスピン1個の上下動を管理するclass Spin を定義する。

//---- declaration and difinition of class Spin
class Spin
{
  int spin;      // the value of spin, it takes +1 or -1
public:
  //-- sets this spin randomly +1 or -1
  void RandSet( void ){
    if( Urand() > 0.5 ) spin = +1; else spin = -1;
  }
  //-- represents the value of spin when this object is called
  operator int (){
    return( spin );
  }
  //-- transits this spin by given probability
  void Transit( int s_sum, double prob ){
    if( Urand() < prob ) spin = +1; else spin = -1;
  }
};
この operator int () を定義すると、このclassのobject sが if( s == +1 ) のように、s に対して int 型の値が 要求されると自動的にs.spin の値が表されることになる。
Transit() methodの引数 s_spin には周りのスピンの総和を 与え、prob にはその状況で上スピンにある確率を与える。

次に格子を管理するclass Lattice を宣言する。

//---- declaration of class Lattice
class Lattice
{
  Spin spin[SIZE][SIZE];    // spins in this lattice
  double prob[5];           // table of transition probability
public:
       Lattice( double T );
  void Sweep( int interval );
  void Measure( Statistic& Energy, Statistic& Magnet );
  void Draw( void );
private:
  inline void Neighbor( int& im, int& ip, int i );
};
先に定義したSpin classのobjectをLattice classが抱える。 小さな配列変数prob[] にはその時の温度での周りのスピンの総和 に対して、上スピンの存在確率を熱浴法により求めた値を設定しておく。 Constructorは、温度を更新した直後に呼び出して prob[] の値の設定と格子の全スピンの向きをランダムに設定させる。
//---- makes transition propability table and let spins randomly turn
Lattice::Lattice( double T )
{
  for( int s_sum=-4 ; s_sum<=+4 ; s_sum+=2 ){
    double eplus = exp( (J*s_sum+H)/T );
    prob[s_sum/2+2] = eplus/(eplus+1.0/eplus);
  }
  for( int ix=0 ; ix<SIZE ; ix++ ){
    for( int iy=0 ; iy<SIZE ; iy++ ){
      spin[ix][iy].RandSet();
    }
  }
}
Sweep() methodは引数 interval 回だけ格子を操引する。 この格子があたかも周期的に並んでいるかのように隣のスピンのindexを 得るための補助method Neighbor も用意しておく。
//---- gets neighbors' indices of given index
inline void Lattice::Neighbor( int& im, int& ip, int i ){
  if( i==0 ) im = SIZE-1; else im = i-1;
  if( i==SIZE-1 ) ip = 0; else ip = i+1;
}

//---- sweeps this lattice [interval] times
void Lattice::Sweep( int interval )
{
  int ix, iy, ixm, ixp, iym, iyp;
  int s_sum;
  
  while( interval-- ){
    for( ix=0 ; ix<SIZE ; ix++ ){
      Neighbor( ixm, ixp, ix );
      for( iy=0 ; iy<SIZE ; iy++ ){
        Neighbor( iym, iyp, iy );
        s_sum = spin[ixp][iy] + spin[ixm][iy] + spin[ix][iyp] + spin[ix][iym];
        spin[ix][iy].Transit( s_sum, prob[s_sum/2+2] );
      }
    }
  }
}
Measure() methodはその時点での格子の全エネルギーと全磁化を 計算して、それを一つの統計量データとして Energy, Magnet のStatistic objectに追加する。
//---- evaluates energy and magnetization of this lattice
void Lattice::Measure( Statistic& Energy, Statistic& Magnet )
{
  int ix, iy, ixm, ixp, iym, iyp;
  int s_sum=0, ss_sum=0;
  
  for( ix=0 ; ix<SIZE ; ix++ ){
    Neighbor( ixm, ixp, ix );
    for( iy=0 ; iy<SIZE ; iy++ ){
      Neighbor( iym, iyp, iy );
      s_sum  += spin[ix][iy];
      ss_sum += spin[ix][iy] * (spin[ixm][iy] + spin[ix][iym]);
    }
  }
  const double energy = -J*ss_sum - H*s_sum;
  const double magnet = (double) s_sum;
  Energy.Add( energy );
  Magnet.Add( magnet );
}
Draw() methodはその時点でのスピンの向きをグラフィックに描く。
//---- draws spins as colored square
void Lattice::Draw( void )
{
  NXClearWindow();

  const int dX = WIN_WIDTH/SIZE, dY = WIN_HEIGHT/SIZE;
  for( int ix=0 ; ix<SIZE ; ix++ ){
    for( int iy=0 ; iy<SIZE ; iy++ ){
      if( spin[ix][iy] == +1 ) NXSetColor( NX_YELLOW ); else NXSetColor( NX_BLUE );
      NXFillRectangle( ix*dX, iy*dY, dX-1, dY-1 );
    }
  }
  NXFlush();
}
さて、計算した結果はグラフィックに描くだけでなく、 生のデータをファイルに保存して、グラフソフトでグラフを 描くことにしよう。4種のデータがあるのでそれぞれにひとつの ファイルを割り当てよう。関数 OpenFiles() は この4個のファイルを開き、その file streamをFILE* 型の 広域配列変数 fptrs[] に格納する。 関数 CloseFiles() はそれらを閉じる。 そして関数 Output() は、ひとつの温度での各量を 画面に出力すると共に各ファイルに書き込む。
//---- opens four files to save data
void OpenFiles( void )
{
  for( int i=0 ; i<4 ; i++ ){
    fptrs[i] = fopen( fnames[i], "wt" );
    if( fptrs[i] == NULL ){
      perror("Error at fopen.");
      exit(1);
    }
  }
}

//---- closes the files
void CloseFiles( void )
{
  for( int i=0 ; i<4 ; i++ ){
    fclose( fptrs[i] );
  }
}

//---- prints and saves each quantities in the four files
void Output( double T, Statistic quantity[] )
{
  printf("\nTemperature    =  %lf\n", T );

  for( int i=0 ; i<4 ; i++ ){
    printf("%-14s = % lf +/- %lf\n", titles[i], quantity[i].Mean(), quantity[i].Error() );
    fprintf( fptrs[i], "%lf  %lf %lf\n", T, quantity[i].Mean(), quantity[i].Error() );
  }
}
指定された温度で実験の総指揮を行う関数 Experiment() である。
//---- executes experiment at given temperature
void Experiment( double T )
{
  Lattice   latt(T);
  Statistic quantity[4];

  latt.Sweep( RELAX );     // makes the lattice equilibrium

  for( int group=0 ; group<GROUPS ; group++ ){
    Statistic Energy0, Magnet0;
    
    // measures fluctuations of energy and magnetization
    for( int sample=0 ; sample<SAMPLES ; sample++ ){
      latt.Sweep( INTERVAL );
      latt.Measure( Energy0, Magnet0 );
    }
    // calculates other data per spin statistically
    quantity[ENERGY].Add( Energy0.Mean()/(SIZE*SIZE) );
    quantity[MAGNET].Add( Magnet0.Mean()/(SIZE*SIZE) );
    quantity[SUSCEP].Add( Magnet0.Variance()/T/(SIZE*SIZE) );
    quantity[SPHEAT].Add( Energy0.Variance()/sqr(T)/(SIZE*SIZE) );
    
    latt.Draw();
  }
  
  Output( T, quantity ); // saves data in files
}
内側のforループでエネルギーと磁化を測定し、その平均値と分散を Energy0, Magnet0 に集めていく。 その分散から帯磁率と比熱を求め、 外側のforループでスピン1個辺りの各種測定量の平均値とエラーを quantity[] に集めていく。

注意すべきはEnergy0, Magnet0 の宣言がループ内にあること である。ループ内で宣言した変数はその回ごとに初期化、解体が行われる ことを思いだそう。つまりプログラムの実行がその宣言文を通るたびに それぞれのconstructorが起動する。このため各回の統計データが 自動的に初期化され独立な実験データとなるのである。

最後にmain() 関数である。

//---- main function
int main(void)
{
  XEvent ev;
  NXOpenWindow( "Ising Model", WIN_WIDTH, WIN_HEIGHT );
  NXEnable( NX_DOUBLEBUFFER );

  OpenFiles();

  for( double T=T_INIT ; T<=T_FINA ; T+=T_STEP ){  // temperature loop
    Experiment( T );
  }

  CloseFiles();
  NXCloseWindow();
  return(0);
}

データをgnuplotでグラフにしたものを下に載せる。 いくつかの温度でのスピンの様子のグラフィックも載せる。


エネルギー


磁化


帯磁率


比熱

相転移温度近くでのスピンの様子



  • 次へ
  • 目次
    Copyright(C) by Naoki Watanabe. Oct 21st, 1995.
    渡辺尚貴 naoki@cms.phys.s.u-tokyo.ac.jp