cleaned up code, added more output info

1 files changed, 48 insertions, 36 deletions

diff --git a/src/new_model.cpp b/src/new_model.cpp
index b96421d..326d530 100644
--- a/src/new_model.cpp
+++ b/src/new_model.cpp
@@ -56,19 +56,21 @@ double vsin(const vector_t<2, double>& v1, const vector_t<2, double>& v2) {
   return v1[0] * v2[1] - v1[1] * v2[0];
 }
 
+// these need to be global so that both wolff hooks have access to them.
+// TODO: investigate using pointers
 double cos_xy;
 vector_t<2, double> sin_xy;
 
 class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i> {
   private:
     double a;
+    double K1;
 
-    unsigned n;
     int A;
 
   public:
     double totalA2;
-    i_measurement(const ising::wolff::system<ising_t, ising_t, graph_i>& S, double a) : a(a) {
+    i_measurement(const ising::wolff::system<ising_t, ising_t, graph_i>& S, double a, double K1) : a(a), K1(K1) {
       A = S.nv;
 
       totalA2 = 0;
@@ -81,7 +83,7 @@ class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i
         A -= s_new - s_old;
       }
 
-      sin_xy[v.prop.sublattice] += a * (s_new - s_old) * vsin(*(v.prop.back), *(v.prop.front));
+      sin_xy[v.prop.sublattice] += K1 * a * (s_new - s_old) * vsin(*(v.prop.back), *(v.prop.front));
 
       cos_xy -= dE;
     }
@@ -93,13 +95,13 @@ class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i
 
 class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, vector_t<2, double>, graph_x> {
   private:
-    double Jxy;
-    double Jxyn;
+    double K1y;
+    double K1yn;
     double a;
     const ising_t* si0;
 
   public:
-    x_measurement(const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>, graph_x>& S, double Jxy, double Jxyn, double a, const ising_t *si0) : Jxy(Jxy), Jxyn(Jxyn), a(a), si0(si0) {
+    x_measurement(const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>, graph_x>& S, double K1y, double K1yn, double a, const ising_t *si0) : K1y(K1y), K1yn(K1yn), a(a), si0(si0) {
     }
 
     void plain_bond_visited(const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>, graph_x>& S, const graph_x::halfedge& e, const vector_t<2, double>& s1_new, double dE) override {
@@ -108,9 +110,9 @@ class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, vec
       vector_t<2, double> dsin = (vsin(s1_new, S.s[e.neighbor.ind]) - vsin(S.s[e.self.ind], S.s[e.neighbor.ind])) * e.prop.dNeighborSelf;
 
       if (e.prop.distance == 0) {
-        sin_xy += (Jxy + (*si0) * a * *(e.prop.spin)) * dsin;
+        sin_xy += K1y * (1 + (*si0) * a * *(e.prop.spin)) * dsin;
       } else if (e.prop.distance == 1) {
-        sin_xy += Jxyn * dsin;
+        sin_xy += K1yn * dsin;
       }
     }
 };
@@ -119,18 +121,19 @@ int main (int argc, char *argv[]) {
 
   // set defaults
   unsigned N = (unsigned)1e4;
-  const unsigned D = 2;
-  unsigned L = 128;
-  double T = 1.0;
-  double Ji = 1.0;
-  double Jx = 1.0;
-  double Jxn = 0.0;
-  double a = 1.5;
+
+  const unsigned D = 2; // we're always in 2D
+  unsigned L = 128;     // system size
+  double T = 1.0;       // temperature
+  double J = 1.0;      // ising nn coupling
+  double K1 = 1.0;      // xy nn coupling
+  double K2 = 0.0;     // xy nnn coupling
+  double α = 1.5;       // ising-xy coupling
 
   int opt;
 
   // take command line arguments
-  while ((opt = getopt(argc, argv, "N:L:T:I:X:a:n:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:L:T:J:K:a:n:")) != -1) {
     switch (opt) {
       case 'N': // number of steps
         N = (unsigned)atof(optarg);
@@ -138,20 +141,20 @@ int main (int argc, char *argv[]) {
       case 'L': // linear size
         L = atoi(optarg);
         break;
-      case 'T': // temperature
+      case 'T':
         T = atof(optarg);
         break;
-      case 'I': // temperature
-        Ji = atof(optarg);
+      case 'J':
+        J = atof(optarg);
         break;
-      case 'X': // temperature
-        Jx = atof(optarg);
+      case 'K':
+        K1 = atof(optarg);
         break;
-      case 'a': // temperature
-        a = atof(optarg);
+      case 'n':
+        K2 = atof(optarg);
         break;
-      case 'n': // temperature
-        Jxn = atof(optarg);
+      case 'a':
+        α = atof(optarg);
         break;
       default:
         exit(EXIT_FAILURE);
@@ -165,16 +168,16 @@ int main (int argc, char *argv[]) {
   // antiferromagnetic Ising coupling
   std::function <double(const ising_t&, const ising_t&)> Zi = [=] (const ising_t& s1, const ising_t s2) -> double {
     if (s1.x == s2.x) {
-      return -Ji;
+      return -J;
     } else {
-      return Ji;
+      return J;
     }
   };
 
   // ising field depends on dot product of surrounding xy spins
   std::function <double(const graph_i::vertex&, const ising_t&)> Bi =
-    [a] (const graph_i::vertex& v, const ising_t& s) -> double {
-      return a * s * (*(v.prop.back) * *(v.prop.front));
+    [K1, α] (const graph_i::vertex& v, const ising_t& s) -> double {
+      return K1 * α * s * (*(v.prop.back) * *(v.prop.front));
     };
 
   // initialize Ising diagonal-lattice graph
@@ -213,11 +216,11 @@ int main (int argc, char *argv[]) {
   // ferromagnetic XY coupling, with nearest-neighbor energy that deponds on
   // state of Ising spin on the bond
   std::function <double(const graph_x::halfedge&, const vector_t<2, double>&, const vector_t<2, double>&)>
-    Zxy = [Jx, Jxn, a, &si] (const graph_x::halfedge& e, const vector_t<2, double>& s1, const vector_t<2, double>& s2) -> double {
+    Zxy = [K1, K2, α, &si] (const graph_x::halfedge& e, const vector_t<2, double>& s1, const vector_t<2, double>& s2) -> double {
       if (e.prop.distance == 0) {
-        return (Jx + si.s0 * a * *(e.prop.spin)) * (s1 * s2);
+        return K1 * (1 + si.s0 * α * *(e.prop.spin)) * (s1 * s2);
       } else if (e.prop.distance == 1) {
-        return Jxn * (s1 * s2);
+        return K1 * K2 * (s1 * s2);
       } else {
         return 0;
       }
@@ -326,11 +329,11 @@ int main (int argc, char *argv[]) {
     si.s[pow(L, D) + i].x = true;
   }
 
-  cos_xy = Jx * 2 * pow(L, 2) + Jxn * 2 * pow(L, 2);
+  cos_xy = K1 * 2 * pow(L, 2) + K2 * 2 * pow(L, 2);
   sin_xy.fill(0);
 
-  i_measurement mi(si, a);
-  x_measurement mxy(sxy, Jx, Jxn, a, &(si.s0));
+  i_measurement mi(si, α, K1);
+  x_measurement mxy(sxy, K1, K2, α, &(si.s0));
 
   double total_cos = 0;
   double total_sin2 = 0;
@@ -344,9 +347,18 @@ int main (int argc, char *argv[]) {
 //    std::cout << cos_xy / sxy.ne << " " << sin_xy / pow(sxy.ne, 2) << " " << E << "\n";
   }
 
+  std::ifstream checkfile("out.dat");
+  bool already_exists = checkfile.good();
+  if (already_exists) {
+    checkfile.close();
+  }
+
   std::ofstream outfile;
   outfile.open("out.dat", std::ios::app);
-  outfile << L << " " << Ji << " " << Jx << " " << a << " " << mi.totalA2 / N / pow(si.nv, 2) << " " << (total_cos - total_sin2 / T) / N / sxy.ne << " " <<(total_cos + total_sin2 / T) / N / sxy.ne << "\n";
+  if (!already_exists) {
+    outfile << "N L J K1 K2 \\[Alpha] M \\[Rho] \\[Sigma]\n";
+  }
+  outfile << N << " " << L << " " << J << " " << K1 << " " << K2 << " " << α << " " << mi.totalA2 / N / pow(si.nv, 2) << " " << (total_cos - total_sin2 / T) / N / sxy.ne << " " <<(total_cos + total_sin2 / T) / N / sxy.ne << "\n";
   outfile.close();
 }