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