Added - algorithms.

diff --git a/.gitignore b/.gitignore
index f6507bdc2b8028b8edb6ca0de030e79baeb0e2e0..9b2b6d79d0e4a27bab67a0b6f11131a16cf825aa 100644 (file)
--- a/.gitignore
+++ b/.gitignore
@@ -8,3 +8,5 @@ heatbath-nor
 heatbath-1mr
 heatbath-r
 heatbath-1mr+
+heatbath-1mr-
+heatbath-r-

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 41c8e074f076cfadb6a35c06ba36212620e72ec8..3e554d45ff163e52c2ec8b8beee5c4e7d6af6d2e 100644 (file)
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -26,3 +26,11 @@ target_link_libraries(heatbath-1mr o815 gsl gslcblas lat_neigh)
 add_executable(heatbath-1mr+ heatbath.cpp)
 set_target_properties(heatbath-1mr+ PROPERTIES COMPILE_DEFINITIONS "ALGORITHM_1MRPLUS")
 target_link_libraries(heatbath-1mr+ o815 gsl gslcblas lat_neigh)
+
+add_executable(heatbath-1mr- heatbath.cpp)
+set_target_properties(heatbath-1mr- PROPERTIES COMPILE_DEFINITIONS "ALGORITHM_1MRMINUS")
+target_link_libraries(heatbath-1mr- o815 gsl gslcblas lat_neigh)
+
+add_executable(heatbath-r- heatbath.cpp)
+set_target_properties(heatbath-r- PROPERTIES COMPILE_DEFINITIONS "ALGORITHM_RMINUS")
+target_link_libraries(heatbath-r- o815 gsl gslcblas lat_neigh)

diff --git a/heatbath.cpp b/heatbath.cpp
index a6f08595d7c0a5a4c22ed655ed19c699166c9f2f..0c1348dbefc1b082b0d19209cfd17cb3ee06b38c 100644 (file)
--- a/heatbath.cpp
+++ b/heatbath.cpp
@@ -14,6 +14,12 @@
 #elif ALGORITHM_1MRPLUS
 #define ALGORITHM "1mr+"
 #include "sim-1mr+.hpp"
+#elif ALGORITHM_1MRMINUS
+#define ALGORITHM "1mr-"
+#include "sim-1mr-.hpp"
+#elif ALGORITHM_RMINUS
+#define ALGORITHM "r-"
+#include "sim-r-.hpp"
 #else
 #error NO UPDATE-ALGORITHM DEFINED!
 #endif

diff --git a/obs_phip2_hist.hpp b/obs_phip2_hist.hpp
index a239d0e444600d1231d5047c85bbd905c98cbfd7..60c86b8e41e40bb18cfbf3152fe7acfdd5c99822 100644 (file)
--- a/obs_phip2_hist.hpp
+++ b/obs_phip2_hist.hpp
@@ -15,7 +15,7 @@ class obs_phip2_hist : public o815::obs {
 
 public:
   struct obsmem {
-    double phip2[2];
+    double phip2[4];
   };
   obs_phip2_hist(o815 *_O815);
   
@@ -28,7 +28,7 @@ private:
 };
 
 obs_phip2_hist::obs_phip2_hist(o815 *_O815) : o815::obs("phip2", 
-                                                       _O815->paraQ->getParaNames() + "nthstep:iruv:phip2absdiffequi", _O815, 
+                                                       _O815->paraQ->getParaNames() + "nthstep:imode:kt:kx:phip2absdiffequi", _O815, 
                                                        sizeof(obsmem), 
                                                        "_hist") {
   OM = (obsmem*)obsMem;
@@ -36,33 +36,30 @@ obs_phip2_hist::obs_phip2_hist(o815 *_O815) : o815::obs("phip2",
 }
 
 void obs_phip2_hist::_meas(bool loadedobs, const int& nthmeas) {
-  if (!loadedobs) {
-    for (int kpiruv=0; kpiruv<2; kpiruv++) {
-      double ppt = 2.*M_PI/O815->comargs.lsize[0] * ( 0 + kpiruv * O815->comargs.lsize[0]/2 ) ;
-      double ppx = 2.*M_PI/O815->comargs.lsize[1] * ( 0 + kpiruv * O815->comargs.lsize[1]/2 );
+  for (int kpimode=0; kpimode<4; kpimode++) {
+    const double ppt = 2.*M_PI/O815->comargs.lsize[0] * ( 0 + int( kpimode/4. * O815->comargs.lsize[0] ) ) ;
+    const double ppx = 2.*M_PI/O815->comargs.lsize[1] * ( 0 + int( kpimode/4. * O815->comargs.lsize[1] ) );
 
-      OM->phip2[ kpiruv ] = 0;
+    if (!loadedobs) {
+      OM->phip2[ kpimode ] = 0;
        
       for (int ixt = 0; ixt < O815->comargs.lsize[0]; ixt++)
        for (int ixx = 0; ixx < O815->comargs.lsize[1]; ixx++)
          for (int ixpt = 0; ixpt < O815->comargs.lsize[0]; ixpt++)
            for (int ixpx = 0; ixpx < O815->comargs.lsize[1]; ixpx++)
-             OM->phip2[ kpiruv ] += real( conj( Sim->conf[ ixt*O815->comargs.lsize[1] + ixx ].phi )
-                                          * Sim->conf[ ixpt*O815->comargs.lsize[1] + ixpx ].phi
-                                          * exp ( _i_*(double)ppx*(double)(ixx-ixpx) + _i_*(double)ppt*(double)(ixt-ixpt) ) );
+             OM->phip2[ kpimode ] += real( conj( Sim->conf[ ixt*O815->comargs.lsize[1] + ixx ].phi )
+                                           * Sim->conf[ ixpt*O815->comargs.lsize[1] + ixpx ].phi
+                                           * exp ( _i_*(double)ppx*(double)(ixx-ixpx) + _i_*(double)ppt*(double)(ixt-ixpt) ) );
        
-      OM->phip2[ kpiruv ] /= Sim->LSIZE2;
+      OM->phip2[ kpimode ] /= Sim->LSIZE2;
     }
-  }
-  
-  for (int kpiruv=0; kpiruv<2; kpiruv++) {
-    double ppt = 2.*M_PI/O815->comargs.lsize[0] * ( 0 + kpiruv * O815->comargs.lsize[0]/2 ) ;
-    double ppx = 2.*M_PI/O815->comargs.lsize[1] * ( 0 + kpiruv * O815->comargs.lsize[1]/2 );
-
+    
     *out << O815->paraQ->getParaVals();
     *out << "\t" << ( Sim->nequi + (nthmeas+1)*Sim->nskip );
-    *out << "\t" << kpiruv;
-    *out << "\t" << abs( OM->phip2[ kpiruv ] - 1./(ppt*ppt + ppx*ppx + Sim->m*Sim->m) ) << endl;
+    *out << "\t" << kpimode;
+    *out << "\t" << 0 + int( kpimode/4. * O815->comargs.lsize[0] );
+    *out << "\t" << 0 + int( kpimode/4. * O815->comargs.lsize[1] );
+    *out << "\t" << abs( OM->phip2[ kpimode ] - 1./(ppt*ppt + ppx*ppx + Sim->m*Sim->m) ) << endl;
   }
 };
 

diff --git a/sim-1mr-.hpp b/sim-1mr-.hpp
new file mode 100644 (file)
index 0000000..00e15be
--- /dev/null
+++ b/sim-1mr-.hpp
@@ -0,0 +1,79 @@
+#ifndef SIM_HPP
+#define SIM_HPP
+
+#include <gsl/gsl_rng.h>
+#include <complex>
+#include <math.h>
+
+#include "latlib/neigh.h"
+
+class sim : public o815::sim {
+public:
+  struct siteconf {
+    complex<double> phi;
+  };
+  sim(o815 *_O815);
+  siteconf* conf;
+  unsigned int LSIZE2;
+  double M;
+  double m;
+  neigh *nb;
+
+private:
+  void _makeSweep();
+  void _newParas();
+  gsl_rng* rangsl;
+  void updatePhi (const int& x);
+};
+
+
+
+sim::sim(o815 *_O815) : o815::sim( _O815, 
+                                  sizeof(siteconf)*
+                                  (_O815->comargs.lsize[0]*_O815->comargs.lsize[1]) ) {
+  conf = (siteconf*)confMem;
+
+  rangsl = gsl_rng_alloc(gsl_rng_ranlxs0);
+  gsl_rng_set(rangsl, time(NULL));
+
+  LSIZE2 = _O815->comargs.lsize[0] * _O815->comargs.lsize[1];
+
+  nb = new neigh(2, _O815->comargs.lsize[0], _O815->comargs.lsize[1]);
+}
+
+void sim::updatePhi (const int& x) 
+{
+  const double theta = gsl_rng_uniform(rangsl) * 2*M_PI;
+  complex<double> V=0;
+
+  for (int nu=0; nu<4; nu++)
+    V += conf[ (*nb)[x*4+nu] ].phi;
+
+  const double r = exp( - M * norm(conf[x].phi)
+                       - 1./M * norm(V)
+                       + 2 * real( conf[x].phi * conj(V) ) );
+
+  double oldarg = arg(conf[x].phi);
+
+  conf[x].phi = sqrt(std::log( 1./(1-r) )) / sqrt(M)
+    * polar(1.0, oldarg)
+    + V / M;
+}
+
+void sim::_makeSweep() {  
+  for (int ichecker=0; ichecker<2; ichecker++)
+    for (int it=0; it<O815->comargs.lsize[0]; it++)
+      for (int iy=(it+ichecker)%2; iy<O815->comargs.lsize[1]; iy+=2)
+       updatePhi( it*O815->comargs.lsize[1] + iy );
+}
+
+void sim::_newParas() {
+  m = (*O815->paraQ)["mass"];
+  M = pow( (*O815->paraQ)["mass"], 2 ) + 4;
+  *log << "SIM: Resetting fields." << endl << flush;
+
+  for (int ix=0; ix<LSIZE2; ix++)
+    conf[ix].phi = 1;
+}
+
+#endif

diff --git a/sim-r-.hpp b/sim-r-.hpp
new file mode 100644 (file)
index 0000000..24f7ca4
--- /dev/null
+++ b/sim-r-.hpp
@@ -0,0 +1,76 @@
+#ifndef SIM_HPP
+#define SIM_HPP
+
+#include <gsl/gsl_rng.h>
+#include <complex>
+#include <math.h>
+
+#include "latlib/neigh.h"
+
+class sim : public o815::sim {
+public:
+  struct siteconf {
+    complex<double> phi;
+  };
+  sim(o815 *_O815);
+  siteconf* conf;
+  unsigned int LSIZE2;
+  double M;
+  double m;
+  neigh *nb;
+
+private:
+  void _makeSweep();
+  void _newParas();
+  gsl_rng* rangsl;
+  void updatePhi (const int& x);
+};
+
+
+
+sim::sim(o815 *_O815) : o815::sim( _O815, 
+                                  sizeof(siteconf)*
+                                  (_O815->comargs.lsize[0]*_O815->comargs.lsize[1]) ) {
+  conf = (siteconf*)confMem;
+
+  rangsl = gsl_rng_alloc(gsl_rng_ranlxs0);
+  gsl_rng_set(rangsl, time(NULL));
+
+  LSIZE2 = _O815->comargs.lsize[0] * _O815->comargs.lsize[1];
+
+  nb = new neigh(2, _O815->comargs.lsize[0], _O815->comargs.lsize[1]);
+}
+
+void sim::updatePhi (const int& x) 
+{
+  const double r = gsl_rng_uniform(rangsl);
+  const double theta = gsl_rng_uniform(rangsl) * 2*M_PI;
+  complex<double> V=0;
+
+  for (int nu=0; nu<4; nu++)
+    V += conf[ (*nb)[x*4+nu] ].phi;
+
+  double oldarg = arg(conf[x].phi);
+
+  conf[x].phi = sqrt(std::log( 1./(1-r) )) / sqrt(M)
+    * polar(1.0, oldarg)
+    + V / M;
+}
+
+void sim::_makeSweep() {  
+  for (int ichecker=0; ichecker<2; ichecker++)
+    for (int it=0; it<O815->comargs.lsize[0]; it++)
+      for (int iy=(it+ichecker)%2; iy<O815->comargs.lsize[1]; iy+=2)
+       updatePhi( it*O815->comargs.lsize[1] + iy );
+}
+
+void sim::_newParas() {
+  m = (*O815->paraQ)["mass"];
+  M = pow( (*O815->paraQ)["mass"], 2 ) + 4;
+  *log << "SIM: Resetting fields." << endl << flush;
+
+  for (int ix=0; ix<LSIZE2; ix++)
+    conf[ix].phi = 1;
+}
+
+#endif