From a46a187b188cdafc6dc22916dd4c8b10fc356af1 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 29 Jan 2019 13:29:12 -0500 Subject: cleaned up code, added more output info --- src/new_model.cpp | 84 +++++++++++++++++++++++++++++++------------------------ 1 file 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 { private: double a; + double K1; - unsigned n; int A; public: double totalA2; - i_measurement(const ising::wolff::system& S, double a) : a(a) { + i_measurement(const ising::wolff::system& 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, 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, 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, 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, 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, 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 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 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 &, 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(); } -- cgit v1.2.3-70-g09d2