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