fixed bugs in Ising, p for bonds pre-computed

1 files changed, 23 insertions, 24 deletions

diff --git a/spheres.cpp b/spheres.cpp
index 106ce4e..b0f4fb4 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -34,20 +34,20 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  std::function<double(spin<double, D, double>, spin<double, D, double>)> Z =
-    [L] (spin<double, D, double> s1, spin<double, D, double> s2) -> double {
-      vector<double, D> diff = s1.x - s2.x;
-      for (unsigned i = 0; i < D; i++) {
-        if (fabs(diff(i)) > L / 2) {
-          diff(i) = L - fabs(diff(i));
-        } else {
-          diff(i) = fabs(diff(i));
-        }
-      }
-      if (diff.transpose() * diff < pow(s1.s + s2.s, 2)) {
-        return -std::numeric_limits<double>::infinity();
+  std::function<double(spin<double, D, double>, spin<double, D, double>, spin<double, D, double>)> Z =
+    [L] (spin<double, D, double> s1, spin<double, D, double> s2, spin<double, D, double> s1_new) -> double {
+      vector<double, D> diff_old = diff(L, s1.x, s2.x);
+      vector<double, D> diff_new = diff(L, s1_new.x, s2.x);
+
+      double rad_sum = pow(s1.s + s2.s, 2);
+
+      bool old_overlap = diff_old.transpose() * diff_old < rad_sum;
+      bool new_overlap = diff_new.transpose() * diff_new < rad_sum;
+
+      if (new_overlap) {
+        return 1.0;
       } else {
-        return 0;
+        return 0.0;
       }
     };
 
@@ -60,16 +60,10 @@ int main(int argc, char* argv[]) {
     [] (model<double, D, double>& m, unsigned i0, spin<double, D, double> s1) -> std::set<unsigned> {
       std::set<unsigned> nn;
       if (i0 < m.s.size()) {
-        std::set<unsigned> os1 = m.dict.on_site(s1.x);
-        std::set<unsigned> nn0 = m.dict.nearest_neighbors(m.s[i0].x);
-        std::set<unsigned> nn1 = m.dict.nearest_neighbors(s1.x);
-        std::set<unsigned> nnn0 = m.dict.next_nearest_neighbors(m.s[i0].x);
-        std::set<unsigned> nnn1 = m.dict.next_nearest_neighbors(s1.x);
-        nn.insert(nn0.begin(), nn0.end());
-        nn.insert(nn1.begin(), nn1.end());
-        nn.insert(nnn0.begin(), nnn0.end());
-        nn.insert(nnn1.begin(), nnn1.end());
-        nn.insert(os1.begin(), os1.end());
+        std::set<unsigned> n_old = m.dict.neighbors(m.s[i0].x, 1);
+        std::set<unsigned> n_new = m.dict.neighbors(s1.x, 1);;
+        nn.insert(n_old.begin(), n_old.end());
+        nn.insert(n_new.begin(), n_new.end());
         nn.insert(m.s.size());
       } else {
         for (unsigned i = 0; i < m.s.size(); i++) {
@@ -94,11 +88,16 @@ int main(int argc, char* argv[]) {
 
   sphere.wolff(T, N, rng);
 
+  std::ofstream snapfile;
+  snapfile.open("sphere_snap.dat");
+
   for (spin<double, D, double> s : sphere.s) {
     spin<double, D, double> rs = sphere.s0.inverse().act(s);
-    std::cout << s.x.transpose() << "\n";
+    snapfile << rs.x.transpose() << "\n";
   }
 
+  snapfile.close();
+
   return 0;
 }