rmc.cpp

   1 /*
   2  * Copyright (c) 2009 Mauro Iazzi
   3  *
   4  * Permission is hereby granted, free of charge, to any person
   5  * obtaining a copy of this software and associated documentation
   6  * files (the "Software"), to deal in the Software without
   7  * restriction, including without limitation the rights to use,
   8  * copy, modify, merge, publish, distribute, sublicense, and/or sell
   9  * copies of the Software, and to permit persons to whom the
  10  * Software is furnished to do so, subject to the following
  11  * conditions:
  12  *
  13  * The above copyright notice and this permission notice shall be
  14  * included in all copies or substantial portions of the Software.
  15  *
  16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  18  * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  20  * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  21  * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  22  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  23  * OTHER DEALINGS IN THE SOFTWARE.
  24  */
  25
  26 #include "qmc.hpp"
  27 #include "config.hpp"
  28
  29 using namespace std;
  30
  31 int main (int argc, char **argv) {
  32     Config conf (argc, argv);
  33     SpinChain init = conf.init;
  34     double alpha = conf.alpha;
  35     int nspins = conf.nspins;
  36     int niter = conf.niter;
  37     std::ofstream &out = conf.out;
  38     bool print_state = conf.print_state;
  39     Energy en(conf.J, conf.D);
  40     int L = conf.L;
  41
  42     int index = 0; // first of reptile
  43     Walker start (VarProb(alpha), init, en);
  44     vector<Walker> reptile(L, start);
  45
  46     double a = 0.0; // acceptance ratio
  47     bool accepted = true; // first value is always accepted
  48     int d = rand()<(RAND_MAX/2)?1:-1;
  49     int first = 0;
  50     int last = L;
  51     double e = 0.0;
  52     double w = 0.0;
  53     std::string state_str = init.toString();
  54     for (int i=0;i<L-1;i++) {
  55         reptile[i].copyOn(&reptile[i+1]);
  56         reptile[i+1].step();
  57     }
  58     for (int i=1;i<=niter;i++) {
  59         d = rand()<(RAND_MAX/2)?1:-1;
  60         first = (L+index+d)%L;
  61         last = (L+index-d)%L;
  62         a = reptile[last].weight()/reptile[first].weight();
  63         accepted = a>(rand()/(RAND_MAX+1.0));
  64         if (!accepted) {
  65             // e still contains the correct value
  66         } else {
  67             //e = en.compute(mc.state(), mc.helper(), mc.probability());
  68             reptile[last].copyOn(&reptile[index]);
  69             //cerr << "accepted d=" << d << " " << reptile[index].step() << endl;
  70             reptile[index].step();
  71             e = ( reptile[index].localEnergy() + reptile[first].localEnergy() ) / 2.0;
  72             if (print_state) state_str = reptile[index].state().toString();
  73             //w = std::exp(start.probability().compute(reptile[(d==1?first+1:index)%L].state()));
  74             //w = w*w;
  75             w = 1.0;
  76             for (int i=0;i<L;i++) {
  77                 w *= (reptile[i].coeff());
  78             }
  79         }
  80         if (conf.M) {
  81                 e = reptile[(L+(index-first)/2)%L].state().staggeredMagnetization();
  82         }
  83         if (out.is_open()) {
  84             if (print_state) {
  85                 out << state_str << " " << (e/nspins) << " " << w << "\n";
  86             } else {
  87                 out << "* " << (e/nspins) << " " << w << "\n";
  88             }
  89         }
  90         if (accepted) index = last;
  91     }
  92     if (out.is_open()) {
  93         out.close();
  94     }
  95     return 0;
  96 }
  97