Tutorials:MonteCarloHOWTO/ja

概要

ALPSアプリケーションでは、すでにいくつかのモンテカルロ（MC）プログラムが入手可能ですが、ALPSの主な目的は、開発者が自身のMCプログラムを簡単に、短時間で開発可能なように手助けをすることです。ALPSアプリケーションには、すべてのMCアルゴリズムで実行されるいくつかの普遍的な機能があり、もう再開発の必要はありません。この章では、ALPSライブラリが提案するシンプルなフレームワークを説明し、古典モンテカルロ、量子モンテカルロ法についていくつかの例で演習をおこないます。

ALPSソースのexample/schedulerディレクトリにイジングモデルのサンプルがあります。このサンプル計算例を用いて、以下で解説します。

　ALPSライブラリを使用する利点　

• 　エラーバーと自己相関時間の自動計算
• 　自動並列化
• 　XML形式による計算結果の自動出力とこれらを扱うコマンドラインツール
• 　負符号に関するリウェイティング（量子モンテカルロシミュレーションの負符号問題）
• 　簡単なチェックポイント作成と、そこからのリスタート機能
• 　入力パラメータの容易な制御
• 　測定の追加が容易
• 　乱数の容易な取扱い
• 　格子の容易な取扱い
• 　量子ハミルトニアンの容易な取扱い（量子モンテカルロ計算）
• ...

わずかなプログラミング/移植によって、これらの利点が得られます。

ALPS MCプログラム使用の第一歩

このチュートリアルは未完成です。コメント、ご提案があれば"comp-phys-alps-users@lists.comp-phys.org"までお願いします。

まず最初に、ALPSライブラリを使用するためにC++言語を使ってください。また、新しいMCアプリケーションの開発をおこなうが、使うべき内部データ構造やアルゴリズムは分かっているものと仮定します。

ALPSライブラリを使うために、ALPSで定義されたクラスから継承した、自分のMC 用のクラスを定義します。今回はこのクラスをMyMonteCarloと呼び、MyMonteCarlo.hppヘッダーファイルに定義します。

 #ifndef MYMC_HPP
 #define MYMC_HPP
 #include <alps/scheduler.h>
 #include <alps/alea.h>
 #include <alps/scheduler/montecarlo.h>
 #include <alps/osiris.h>
 #include <alps/osiris/dump.h>
 #include <alps/expression.h>
 using namespace std;
 class MyMonteCarlo : public alps::scheduler::MCRun{
 public :
   MyMonteCarlo(const alps::ProcessList &,const alps::Parameters &,int);
   static void print_copyright(std::ostream &);
   void save(alps::ODump &) const;
   void load(alps::IDump &);
   void dostep();
   bool is_thermalized() const;
   double work_done() const;
 private :
 // your own internal data here ...
 };
 #endif

各行が何を意味するか説明していきます。

• 　最初に、ALPSのヘッダーファイルをインクルードする必要があります。（他のALPS関数を利用する場合は、さらに追加する必要あり）
• 　ALPSで定義されているクラスから派生させることによってMyMonteCarloクラスを定義します。上記のヘッダーファイルでは、最も簡単なクラスであるalps::csheduler::MCRunクラスを使用することを意味しています。もしALPSの格子機能を使用したい場合は、次のように置き換えてください。

 class MyMonteCarlo : public alps::scheduler::LatticeMCRun<>{

ヘッダーの追加もおこなってください。

 #include <alps/lattice.h>

もし、さらにモデルライブラリ（格子量子モデル）を使用する場合は、次の様に置き換えてください。

 class MyMonteCarlo : public alps::scheduler::LatticeModelMCRun<>{

• 　publicにある、MyMonteCarlo(const alps::ProcessList&,const alps::Parameters&,int)はコンストラクタです。そして、すべてのパラメータの初期化が必要でしょう。
• 　print_copyright(std::ostream&)は、各シミュレーションの最初に、任意の情報を表示する関数です。
• 　save,loadは、各チェックポイントの後で再スタートするために、何がディスク上にsave(load)する必要があるか決定する関数です。
• 　dostep関数はMCの中心となる関数です。毎回のMCステップごとに実行されます。
• 　is_thermalizedは熱平衡化が終わっているかどうかをALPSライブラリに伝える関数です。
• 　work_doneは、シミュレーションの何パーセントが実行済みかALPSライブラリに伝える関数です。　

最後に、これらの関数を正しく呼びだすプログラムを、MyMonteCarlo.cppとして書きます。

ALPS MCプログラムの作成

copyright関数

まず、MyMonteCarlo.cppの中で最も簡単なprint_copyright()関数から始めます。

 #include "MyMonteCarlo.hpp";
 /************************************ ALPS functions **********************************************/
 // Copyright statement
 void MyMonteCarlo::print_copyright(std::ostream & out)
 {
   out << "My own ALPS Monte Carlo program v. 0.1\n"
       << "  copyright (c) 2006 by Myself,\n"
       << "  available from the author on request\n\n";
 }

モンテカルロステップの制御

次に、MCステップについて説明します。典型的な例は次の通りです。まず、平衡化のためにある一定のステップ数だけ更新し、続いて、あらかじめ決めておいた回数だけ更新しながら物理量を測定します。普通、各測定の間に何回かのモンテカルロステップを実行します。これは、内部データ構造で、次の変数(in MyMonteCarlo.hpp)を設定することで実現できます。

 private :
  // your own internal data here ...
  int Nb_Steps; 
  int Nb_Thermalisation_Steps; 
  int Each_Measurement;
  int Steps_Done_Total; 
  int Measurements_Done;
  void do_update();
  void do_measurements();
  // the rest of your own internal data here ...
 };

Nb_Thermalisation_Stepsは平衡化ステップ数、Nb_Stepsは平衡化後のMCステップ数です。Each_Measuremenは各測定のステップ数をあらわします。これらの3つの数字はコンストラクタで初期化します。Steps_Done_Totalは、今までに平衡化も含めて何回更新したのかを保存します。最後に、Measurements_Doneは前回の測定から何回更新したのかを保存します。do_update()とdo_measurements()関数はそれぞれ毎回のMC更新時および測定時に実行されます。

 bool MyMonteCarlo::is_thermalized() const
   {  return (Steps_Done_Total >= Nb_Thermalisation_Steps); }

この関数は、現在のMCステップ数が平衡化ステップ数より大きいとき、trueを返し、それ以外はfalseを返します。

 double MyMonteCarlo::work_done() const
 { return (is_thermalized() ? (Steps_Done_Total-Nb_Thermalisation_Steps)/double(Nb_Steps) :0.); }

シミュレーションが平衡化していなければ0を返し、そうでなければ測定ステップのうちどれだけが実行されたかを返します。最後に、dostep()関数は次のようになります。

 void MyMonteCarlo::dostep()
 { do_update(); // you'll have to define what this function does later
   ++Steps_Done_Total; // increment the number of steps done
   if (is_thermalized()) // do measurements only if simulation thermalized
     { if (++Measurements_Done==Each_Measurement) // do a measurement every Each_Measurement
       { Measurements_Done=0;
       do_measurements(); // you'll have to define what this function does later
       }
    }
 }

もちろん、あなたのMC計算の、1回のMCステップの間に行われることは、あなた自身で定義しなければなりません。それをdo_update() 関数の中に実装してください。 測定は、後述するdo_measurement()関数の中でまとめて行います。

セーブ、ロード関数

ALPSには計算中の内部データをディスクにチェックポイントとして生成する機能があります。今回は、チェックポイントを利用して計算を再スタートするために、整数変数と倍精度浮動小数点数の配列が必要になると仮定します。

 private :
  // the rest of your own internal data here ...
  int Number_of_Spins; 
  int MyOwnVariable;
  std::vector<double> SpinArray;
  ...
 };

save関数は簡単にかけて、

 void MyMonteCarlo::save(alps::ODump& dump) const
   { dump <<  Number_of_Spins << MyOwnVariable << SpinArray;}

load関数も対となるsaveから簡単に実装できます。

 void MyMonteCarlo::load(alps::IDump& dump)
   { dump >>  Number_of_Spins >> MyOwnVariable >> SpinArray;}

ここで、ダンプの方法を説明します。一般的な型（int,double,boolなど）や一般的なコンテナ（vectorやset）が対応可能です。

次のような構造体を定義したとします。

 struct Vertex
   { int vertex_type;
    std::vector<double> coordinates;
    int SomeOtherVariable; }

この構造体のインスタンスを、自分のモンテカルロクラスの使って保存したいと考えると、

 private :
    // the rest of your own internal data here ...
    Vertex MyVertex; 
    ...
   };

構造体の中で、何を保存、読み出しするのかをALPSに指示しなくてはいけません。

 alps::ODump& operator<<(alps::ODump& dump, const Vertex& v)
 { return dump << v.vertex_type << v.coordinates; }

 alps::IDump& operator>>(alps::IDump& dump, Vertex& v)
 { return dump >> v.vertex_type >> v.coordinates;}

こうすることで、ALPSのsave,load関数にVertex型の変数MyVertexを追加することが容易にできます。

 void MyMonteCarlo::save(alps::ODump& dump) const
   { dump <<  Number_of_Spins << MyOwnVariable << SpinArray << MyVertex;}

もし、上記のようにモンテカルロステップ数を設定したならば、save/load関数は次の用に加えることができます。

 void MyMonteCarlo::save(alps::ODump& dump) const
   { dump <<  Number_of_Spins << MyOwnVariable << SpinArray << MyVertex;
     dump <<  Nb_Steps << Measurements_Done; }

do_measurements()関数

この関数では、物理量を測定し、ALPSに引渡します。

 void MyMonteCarlo::do_measurements()
   { double Energy; std::valarray<double> Correl(L);
     // do the measurements in your code (update the Energy and Correl variable) ...
     ...

     // give them to ALPS
     measurements["Energy"] << Energy;
     measurements["Spin Correlations"] << Correl;
   }

わかりますか？驚いたかもしれませんね。ALPSは "Energy"や"Spin Correlations"をどうやって知るのか？また、スカラー値（エネルギーなど）とベクター値（相関関数など）をどうやって区別するのか？この２つの疑問は、クラスのコンストラクタで以下の用に対処されます。

他では、なぜベクター計測値にstd::vectorではなく、std::valarrayを使用しているのか？これはALPS内部の問題に起因しています。ALPSは各測定値が常に同じサイズであることを要請します。