u1casc-ordinary/sim.hpp

   1 #ifndef SIM_HPP
   2 #define SIM_HPP
   3
   4 #include <gsl/gsl_rng.h>
   5 #include <complex>
   6 #include <math.h>
   7
   8 #include "latlib/neigh.h"
   9
  10 #define EPSILONU 1
  11 #define EPSILONPHI 0.5
  12
  13 class sim : public o815::sim {
  14 public:
  15   sim(o815 *_O815);
  16   unsigned int lsize4;
  17   neigh *nb;
  18   complex<double> *U, *phi;
  19   double kappa[2], lambda[2], beta;
  20
  21 private:
  22   void _makeSweep();
  23   void _newParas();
  24   gsl_rng* rangsl;
  25   double rhoPhi(const int& iphi, const int& x0, const complex<double>& candPhi);
  26   double rhoU(const int& x0, const int& nu0, const complex<double>& candU);
  27   int updatePhi(const int& iphi, const int& x0);
  28   int updateU(const int& x0, const int& nu0);
  29 };
  30
  31 sim::sim(o815 *_O815) : o815::sim( _O815,
  32                                    sizeof(complex<double>)*
  33                                    _O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[1]*(2+4) ) {
  34
  35   lsize4 = _O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[1];
  36
  37   nb = new neigh(4, _O815->comargs.lsize[0], _O815->comargs.lsize[0], _O815->comargs.lsize[0], _O815->comargs.lsize[1]);
  38
  39   phi = (complex<double>*)confMem;
  40   U = (complex<double>*)(confMem + sizeof(complex<double>)*lsize4*2);
  41
  42   rangsl = gsl_rng_alloc(gsl_rng_ranlxs0);
  43   gsl_rng_set(rangsl, time(NULL));
  44 }
  45
  46 void sim::_makeSweep() {
  47   for( int ix=0; ix<lsize4; ix++ ) {
  48     for( int inu=0; inu<4; inu++) updateU(ix, inu);
  49     for( int iphi=0; iphi<2; iphi++) updatePhi(iphi, ix);
  50   }
  51 }
  52
  53 void sim::_newParas() {
  54   kappa[0] = (*O815->paraQ)["kappaone"];
  55   kappa[1] = (*O815->paraQ)["kappatwo"];
  56   lambda[0] = (*O815->paraQ)["lambdaone"];
  57   lambda[1] = (*O815->paraQ)["lambdatwo"];
  58   beta = (*O815->paraQ)["beta"];
  59
  60   for(int ix=0; ix<lsize4; ix++) {
  61     for(int i=0; i<2; i++) phi[ i*lsize4 + ix ] = 0;
  62     for(int nu=0; nu<4; nu++) U[ ix*4 + nu ] = 1;
  63   }
  64 }
  65
  66 int sim::updateU(const int& x0, const int& nu0)
  67 {
  68   complex<double> candU = U[x0*4+nu0] * polar(1.0, 2*EPSILONU*( 0.5 - gsl_rng_uniform(rangsl) ));
  69
  70   if ( gsl_rng_uniform(rangsl) <= rhoU(x0, nu0, candU) ) {
  71     U[x0*4 + nu0] = candU;
  72     return 1;
  73   }
  74
  75   return 0;
  76 }
  77
  78 int sim::updatePhi(const int& iphi, const int& x0)
  79 {
  80   complex<double> candPhi = phi[ iphi*lsize4 + x0 ] +
  81     complex<double>(2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) ),
  82                     2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) ));
  83
  84   if ( gsl_rng_uniform(rangsl) <= rhoPhi(iphi, x0, candPhi) ) {
  85     phi[ iphi*lsize4 + x0 ] = candPhi;
  86     return 1;
  87   }
  88
  89   return 0;
  90 }
  91
  92 double sim::rhoPhi(const int& iphi, const int& x0, const complex<double>& candPhi)
  93 {
  94   double deltaS=0;
  95
  96   for( int mu=0; mu<4; mu++) {
  97     if( iphi == 0 ) {
  98       deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * U[ x0*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
  99       deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * conj(U[ (*nb)[x0*8+mu+4]*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
 100       deltaS -= 2 * real( conj(candPhi) * U[ x0*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
 101       deltaS -= 2 * real( conj(candPhi) * conj(U[ (*nb)[x0*8+mu+4]*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
 102     }
 103     else if( iphi == 1 ) {
 104       deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * conj(U[ x0*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
 105       deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * U[ (*nb)[x0*8+mu+4]*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
 106       deltaS -= 2 * real( conj(candPhi) * conj(U[ x0*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
 107       deltaS -= 2 * real( conj(candPhi) * U[ (*nb)[x0*8+mu+4]*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
 108     }
 109   }
 110
 111   deltaS -= kappa[iphi] * norm(phi[ iphi*lsize4 + x0 ]);
 112   deltaS += kappa[iphi] * norm(candPhi);
 113
 114   deltaS -= lambda[iphi] * pow(norm(phi[ iphi*lsize4 + x0 ]),2);
 115   deltaS += lambda[iphi] * pow(norm(candPhi),2);
 116
 117   return exp(-deltaS);
 118 }
 119
 120 double sim::rhoU(const int& x0, const int& nu0, const complex<double>& candU)
 121 {
 122   double deltaS=0;
 123
 124   for( int nu=0; nu<4; nu++ ) {
 125     if( nu == nu0 ) continue;
 126     deltaS += beta * real( U[x0*4+nu0] * U[ (*nb)[x0*8+nu0]*4 + nu ] * conj(U[ (*nb)[x0*8+nu]*4 + nu0 ]) * conj(U[ x0*4 + nu ]) );
 127     deltaS += beta * real( U[ (*nb)[x0*8+nu+4]*4 + nu0 ] * U[ (*nb)[ (*nb)[x0*8+nu+4]*8+nu0 ]*4 + nu ] * conj(U[ x0*4 + nu0 ]) * conj(U[ (*nb)[x0*8+nu+4]*4 + nu ]) );
 128     deltaS -= beta * real( candU * U[ (*nb)[x0*8+nu0]*4 + nu ] * conj(U[ (*nb)[x0*8+nu]*4 + nu0 ]) * conj(U[ x0*4 + nu ]) );
 129     deltaS -= beta * real( U[ (*nb)[x0*8+nu+4]*4 + nu0 ] * U[ (*nb)[ (*nb)[x0*8+nu+4]*8+nu0 ]*4 + nu ] * conj(candU) * conj(U[ (*nb)[x0*8+nu+4]*4 + nu ]) );
 130   }
 131
 132   deltaS += 2 * real( conj(phi[ 0*lsize4 + x0 ]) * U[ x0*4 + nu0 ] * phi[ 0*lsize4 + (*nb)[x0*8+nu0] ]  );
 133   deltaS -= 2 * real( conj(phi[ 0*lsize4 + x0 ]) * candU * phi[ 0*lsize4 + (*nb)[x0*8+nu0] ]  );
 134
 135   deltaS += 2 * real( conj(phi[ 1*lsize4 + x0 ]) * conj(U[ x0*4 + nu0 ]) * phi[ 1*lsize4 + (*nb)[x0*8+nu0] ]  );
 136   deltaS -= 2 * real( conj(phi[ 1*lsize4 + x0 ]) * conj(candU) * phi[ 1*lsize4 + (*nb)[x0*8+nu0] ]  );
 137
 138   return exp(-deltaS);
 139 }
 140
 141 #endif