added examples for no field and site depedence

4 files changed, 201 insertions, 2 deletions

diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 634c600..fcb6c0f 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -2,6 +2,8 @@
 add_executable(ising ising.cpp)
 add_executable(ising_animation ising_animation.cpp)
 add_executable(ising_standalone ising_standalone.cpp)
+add_executable(ising_random_field ising_random_field.cpp)
+add_executable(ising_no_field ising_no_field.cpp)
 add_executable(xy On.cpp)
 add_executable(potts_3 potts.cpp)
 add_executable(clock_5 clock.cpp)
@@ -15,6 +17,8 @@ target_compile_definitions(clock_5 PUBLIC WOLFF_POTTSQ=5)
 target_link_libraries(ising wolff)
 target_link_libraries(ising_animation wolff GL GLU glut)
 target_link_libraries(ising_standalone wolff)
+target_link_libraries(ising_random_field wolff)
+target_link_libraries(ising_no_field wolff)
 target_link_libraries(xy wolff)
 target_link_libraries(potts_3 wolff)
 target_link_libraries(clock_5 wolff)
diff --git a/examples/ising_no_field.cpp b/examples/ising_no_field.cpp
new file mode 100644
index 0000000..eebd672
--- /dev/null
+++ b/examples/ising_no_field.cpp
@@ -0,0 +1,78 @@
+
+#include <getopt.h>
+#include <iostream>
+#include <chrono>
+
+#define WOLFF_NO_FIELD
+#include "simple_measurement.hpp"
+
+#include <wolff/models/ising.hpp>
+#include <wolff/finite_states.hpp>
+
+#include <wolff.hpp>
+
+using namespace wolff;
+
+int main(int argc, char *argv[]) {
+
+  // set defaults
+  N_t N = (N_t)1e4;
+  D_t D = 2;
+  L_t L = 128;
+  double T = 2.26918531421;
+
+  int opt;
+
+  // take command line arguments
+  while ((opt = getopt(argc, argv, "N:D:L:T:")) != -1) {
+    switch (opt) {
+    case 'N': // number of steps
+      N = (N_t)atof(optarg);
+      break;
+    case 'D': // dimension
+      D = atoi(optarg);
+      break;
+    case 'L': // linear size
+      L = atoi(optarg);
+      break;
+    case 'T': // temperature
+      T = atof(optarg);
+      break;
+    default:
+      exit(EXIT_FAILURE);
+    }
+  }
+
+  // define the spin-spin coupling
+  std::function <double(const ising_t&, const ising_t&)> Z = [] (const ising_t& s1, const ising_t& s2) -> double {
+    return (double)(s1 * s2);
+  };
+
+  // initialize the lattice
+  graph G(D, L);
+
+  // initialize the system
+  system<ising_t, ising_t> S(G, T, Z);
+
+  // initialize the random number generator
+  auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
+  std::mt19937 rng{seed};
+
+  // define function that generates self-inverse rotations
+  std::function <ising_t(std::mt19937&, const system<ising_t, ising_t>&, v_t)> gen_r = gen_ising;
+
+  // initailze the measurement object
+  simple_measurement A(S);
+
+  // run wolff N times
+  S.run_wolff(N, gen_r, A, rng);
+
+  // print the result of our measurements
+  std::cout << "Wolff complete!\nThe average energy per site was " << A.avgE() / S.nv
+    << ".\nThe average magnetization per site was " << A.avgM() / S.nv
+    << ".\nThe average cluster size per site was " << A.avgC() / S.nv << ".\n";
+
+  // exit
+  return 0;
+}
+
diff --git a/examples/ising_random_field.cpp b/examples/ising_random_field.cpp
new file mode 100644
index 0000000..b37b0ce
--- /dev/null
+++ b/examples/ising_random_field.cpp
@@ -0,0 +1,93 @@
+
+#include <getopt.h>
+#include <iostream>
+#include <chrono>
+
+#define WOLFF_SITE_DEPENDENCE
+
+#include "simple_measurement.hpp"
+
+#include <wolff/models/ising.hpp>
+
+#include <wolff.hpp>
+
+using namespace wolff;
+
+int main(int argc, char *argv[]) {
+
+  // set defaults
+  N_t N = (N_t)1e4;
+  D_t D = 2;
+  L_t L = 128;
+  double T = 2.26918531421;
+  double H = 0.0;
+
+  int opt;
+
+  // take command line arguments
+  while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
+    switch (opt) {
+    case 'N': // number of steps
+      N = (N_t)atof(optarg);
+      break;
+    case 'D': // dimension
+      D = atoi(optarg);
+      break;
+    case 'L': // linear size
+      L = atoi(optarg);
+      break;
+    case 'T': // temperature
+      T = atof(optarg);
+      break;
+    case 'H': // external field
+      H = atof(optarg);
+      break;
+    default:
+      exit(EXIT_FAILURE);
+    }
+  }
+
+  // define the spin-spin coupling
+  std::function <double(const ising_t&, const ising_t&)> Z = [] (const ising_t& s1, const ising_t& s2) -> double {
+    return (double)(s1 * s2);
+  };
+
+  // initialize the lattice
+  graph G(D, L);
+
+  // initialize the random number generator
+  auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
+  std::mt19937 rng{seed};
+
+  // define the spin-field coupling
+  std::vector<double> H_vals(G.nv);
+  std::normal_distribution<double> distribution(0.0, H);
+  for (v_t i = 0; i < G.nv; i++) {
+    H_vals[i] = distribution(rng);
+  }
+
+  std::function <double(v_t, const ising_t&)> B = [&] (v_t i, const ising_t& s) -> double {
+    return H_vals[i] * s;
+  };
+
+  // initialize the system
+  system<ising_t, ising_t> S(G, T, Z, B);
+
+  // define function that generates self-inverse rotations
+  std::function <ising_t(std::mt19937&, const system<ising_t, ising_t>&, v_t)> gen_r = gen_ising;
+
+  // initailze the measurement object
+  simple_measurement A(S);
+
+  // run wolff N times
+  S.run_wolff(N, gen_r, A, rng);
+
+  // print the result of our measurements
+  std::cout << "Wolff complete!\nThe average energy per site was " << A.avgE() / S.nv
+    << ".\nThe average magnetization per site was " << A.avgM() / S.nv
+    << ".\nThe average cluster size per site was " << A.avgC() / S.nv << ".\n";
+
+  // exit
+  return 0;
+}
+
diff --git a/examples/simple_measurement.hpp b/examples/simple_measurement.hpp
index 518631c..f98cfe6 100644
--- a/examples/simple_measurement.hpp
+++ b/examples/simple_measurement.hpp
@@ -20,7 +20,27 @@ class simple_measurement : public measurement<R_t, X_t> {
     simple_measurement(const system<R_t, X_t>& S) {
       n = 0;
       M = S.nv * S.s[0];
-      E = - (S.ne * S.Z(S.s[0], S.s[0]) + S.nv * S.B(S.s[0]));
+      E = 0;
+
+#ifdef WOLFF_BOND_DEPENDENCE
+      for (v_t i = 0; i < S.nv; i++) {
+        for (v_t j : S.G.adjacency[i]) {
+          E -= 0.5 * S.Z(i, S.s[i], j, S.s[j]);
+        }
+      }
+#else
+      E -= S.ne * S.Z(S.s[0], S.s[0]);
+#endif
+
+#ifndef WOLFF_NO_FIELD
+#ifdef WOLFF_SITE_DEPENDENCE
+      for (v_t i = 0; i < S.nv; i++) {
+        E -= S.B(i, S.s[i]);
+      }
+#else
+      E -= S.nv * S.B(S.s[0]);
+#endif
+#endif
 
       totalE = 0;
       totalM = 0.0 * (S.s[0]);
@@ -40,8 +60,12 @@ class simple_measurement : public measurement<R_t, X_t> {
       M += s_new - s_old;
     }
 
-    void plain_site_transformed(const system<R_t, X_t>&, v_t,  const X_t&) {
+    void plain_site_transformed(const system<R_t, X_t>& S, v_t i, const X_t& si_new) {
       C++;
+
+#ifdef WOLFF_NO_FIELD
+      M += si_new - S.s[i];
+#endif
     }
 
     void ghost_site_transformed(const system<R_t, X_t>&, const R_t&) {}