From b96c82a1c1d95cfa090c4c99e4b4851a5c96e202 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias <jaron@kent-dobias.com> Date: Tue, 25 Feb 2020 16:49:40 -0500 Subject: Started implementing dimers on a torus --- dimers_torus.cpp | 137 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ spheres.cpp | 109 +------------------------------------------ spheres.hpp | 19 ++++++++ 3 files changed, 158 insertions(+), 107 deletions(-) create mode 100644 dimers_torus.cpp diff --git a/dimers_torus.cpp b/dimers_torus.cpp new file mode 100644 index 0000000..f304564 --- /dev/null +++ b/dimers_torus.cpp @@ -0,0 +1,137 @@ + +#include <chrono> +#include <fstream> +#include <iostream> + +#include "animation.hpp" +#include "space_wolff.hpp" +#include "spheres.hpp" +#include "torus_symmetries.hpp" + +const unsigned D = 2; +typedef Model<double, D, TorusGroup<double, D>, Dimer<double, 2>> model; + +Gen<double, D, TorusGroup<double, D>, Dimer<double, D>> eGen(double L) { + std::vector<Vector<double, 2>> torusVectors = torus_vecs<double, 2>(L); + std::vector<Matrix<double, 2>> torusMatrices = torus_mats<double, 2>(); + return [L, torusVectors, + torusMatrices](Model<double, D, TorusGroup<double, D>, Dimer<double, D>>& M, + Rng& r) -> Transformation<double, D, TorusGroup<double, D>, Dimer<double, D>>* { + Matrix<double, 2> m; + Vector<double, 2> t; + + m = r.pick(torusMatrices); + t(0) = r.uniform<double>(0, L); + t(1) = r.uniform<double>(0, L); + t = t - m * t; + + TorusGroup<double, 2> g = TorusGroup<double, 2>({(double)L, t, m}); + + Spin<double, 2, Dimer<double, 2>>* ss = r.pick(M.s); + + return new SpinFlip<double, 2, TorusGroup<double, 2>, Dimer<double, 2>>(M, g, ss); + }; +} + +int main(int argc, char* argv[]) { + const unsigned D = 2; + + double L = 32; + unsigned N = 1000; + double T = 2.0 / log(1.0 + sqrt(2.0)); + double H = 1.0; + unsigned n = 25; + unsigned wait = 1000; + + double k = 1e2; + double a = 0.1; + + int opt; + + while ((opt = getopt(argc, argv, "n:N:L:T:H:a:k:w:")) != -1) { + switch (opt) { + case 'n': + n = (unsigned)atof(optarg); + break; + case 'N': + N = (unsigned)atof(optarg); + break; + case 'L': + L = atof(optarg); + break; + case 'T': + T = atof(optarg); + break; + case 'H': + H = atof(optarg); + break; + case 'a': + a = atof(optarg); + break; + case 'k': + k = atof(optarg); + break; + case 'w': + wait = atoi(optarg); + break; + default: + exit(1); + } + } + + auto zSingle = zSpheresTorus<D>(L, a, k); + + std::function<double(const Spin<double, D, Dimer<double, D>>&, + const Spin<double, D, Dimer<double, D>>&)> + Z = [L, zSingle, a, k](const Spin<double, D, Dimer<double, D>>& s1, + const Spin<double, D, Dimer<double, D>>& s2) -> double { + Spin<double, D, Radius> s11 = {.x = s1.x + s1.s.relativePosition, .s = s1.s.radius}; + Spin<double, D, Radius> s12 = {.x = s1.x - s1.s.relativePosition, .s = s1.s.radius}; + Spin<double, D, Radius> s21 = {.x = s2.x + s2.s.relativePosition, .s = s2.s.radius}; + Spin<double, D, Radius> s22 = {.x = s2.x - s2.s.relativePosition, .s = s2.s.radius}; + + return zSingle(s11, s21) + zSingle(s12, s21) + zSingle(s11, s22) + zSingle(s12, s22); + }; + + std::function<double(Spin<double, D, Dimer<double, D>>)> B = [L, H](Spin<double, D, Dimer<double, D>> s) -> double { + return H * s.x(1); + }; + + auto g1 = eGen(L); + + auto tag = std::chrono::high_resolution_clock::now(); + + Animation<double, D, TorusGroup<double, D>, Dimer<double, D>> A(L, 750, argc, argv, wait, true); + model sphere(L, Z, B); + + Rng rng; + + sphere.s.resize(n); + + unsigned nx = floor(sqrt(n)); + for (unsigned i = 0; i < sphere.s.size(); i++) { + Spin<double, 2, Dimer<double, D>>* ss = new Spin<double, 2, Dimer<double, D>>(); + ss->x = {(i / nx) * L / nx - L / 2, (i % nx) * L / nx - L / 2}; + ss->s = Dimer<double, D>{.relativePosition = {0.2, 0}, .radius = 0.25}; + sphere.s[i] = ss; + sphere.dict.insert(ss); + } + + sphere.wolff(T, {g1}, A, N); + + /* + 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); + snapfile << rs.s << " " << rs.x.transpose() << "\n"; + delete s; + } + + snapfile.close(); + */ + + return 0; +} + diff --git a/spheres.cpp b/spheres.cpp index 4235d66..0067cc8 100644 --- a/spheres.cpp +++ b/spheres.cpp @@ -5,116 +5,11 @@ #include "space_wolff.hpp" #include "torus_symmetries.hpp" +#include "animation.hpp" const unsigned D = 2; typedef Model<double, D, TorusGroup<double, D>, double> model; -class animation : public measurement<double, 2, TorusGroup<double, D>, double> { -private: - uint64_t t1; - uint64_t t2; - unsigned n; - unsigned tmp; - unsigned wait; - double L; - TorusGroup<double, D> s0_tmp; - std::ofstream& file; - -public: - animation(double L, unsigned w, int argc, char* argv[], unsigned wait, std::ofstream& file) : s0_tmp(0), wait(wait), L(50 * L), file(file) { - t1 = 0; - t2 = 0; - n = 0; - - glutInit(&argc, argv); - glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); - glutInitWindowSize(w, w); - glutCreateWindow("wolffWindow"); - glClearColor(0.0, 0.0, 0.0, 0.0); - glMatrixMode(GL_PROJECTION); - glLoadIdentity(); - gluOrtho2D(-1, L+1, -1, L+1); - } - - void pre_cluster(const model& m, unsigned, const Transformation<double, D, TorusGroup<double, D>, double>* t) override { - s0_tmp = m.s0.inverse(); - tmp = 0; - glClearColor(1.0f, 1.0f, 1.0f, 1.0f); - glClear(GL_COLOR_BUFFER_BIT); - glColor3d(1.0, 0.0, 0.0); - glBegin(GL_LINES); - Vector<double, 2> r = t->r.t / 2; - double θ = atan2(r(1), r(0)); - Vector<double, 2> v1 = s0_tmp.act({r(0) + L * sin(θ), r(1) - L * cos(θ)}); - Vector<double, 2> v2 = s0_tmp.act({r(0) - L * sin(θ), r(1) + L * cos(θ)}); - glVertex2d(v1(0), v1(1)); - glVertex2d(v2(0), v2(1)); - if (n > wait) - file << v1.transpose() << " " << v2.transpose() << std::endl; - glEnd(); - for (const Spin<double, 2, double>* s : m.s) { - glBegin(GL_POLYGON); - unsigned n_points = 50; - glColor3f(0.0f, 0.0f, 0.0f); - if (n > wait) - file << s << " " << s->s << " " << s0_tmp.act(*s).x.transpose() << " "; - for (unsigned i = 0; i < n_points; i++) { - glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), - s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); - } - glEnd(); - } - if (n > wait) - file << std::endl; - glLineWidth(3); - glFlush(); - } - - void plain_site_transformed(const model& m, const Transformation<double, D, TorusGroup<double, D>, double>& t) override { - if (n > wait) { - if (t.current().size() > 1) { - glColor3f(0.0f, 0.0f, 1.0f); - } else { - glColor3f(0.0f, 1.0f, 0.0f); - } - for (const Spin<double, 2, double>* s : t.current()) { - if (s != NULL) { - file << s << " " << s0_tmp.act(t.r.act(*s)).x.transpose() << " "; - glBegin(GL_POLYGON); - unsigned n_points = 50; - for (unsigned i = 0; i < n_points; i++) { - glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), - s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); - } - glEnd(); - } else { - file << s << " " << s0_tmp.act({0, 0}).transpose() << " "; - } - } - } - tmp++; - } - - void post_cluster(const model& m) override { - t1 += tmp; - t2 += tmp * tmp; - if (n > wait) { - file << std::endl; - glFlush(); - sleep(2); - } - n++; - } - - void clear() { - t1 = 0; - t2 = 0; - n = 0; - } - - double var() { return (t2 - t1 * t1 / (double)n) / (double)n; } -}; - Gen<double, D, TorusGroup<double, D>, double> eGen(double L) { std::vector<Vector<double, 2>> torusVectors = torus_vecs<double, 2>(L); std::vector<Matrix<double, 2>> torusMatrices = torus_mats<double, 2>(); @@ -196,7 +91,7 @@ int main(int argc, char* argv[]) { auto g = eGen(L); std::ofstream ofile("test.dat"); - animation A(L, 750, argc, argv, 1000, ofile); + Animation<double, D, TorusGroup<double, D>, Radius> A(L, 750, argc, argv, 1000, true); model sphere(L, Z, B); randutils::mt19937_rng rng; diff --git a/spheres.hpp b/spheres.hpp index 040313c..3dc4a5e 100644 --- a/spheres.hpp +++ b/spheres.hpp @@ -20,6 +20,25 @@ zSpheres(double a, double k) { }; } +template <int D> +std::function<double(const Spin<double, D, Radius>&, const Spin<double, D, Radius>&)> +zSpheresTorus(double L, double a, double k) { + return [L, a, k](const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double { + Vector<double, D> d = diff(L, s1.x, s2.x); + + double σ = s1.s + s2.s; + double δ = σ - sqrt(d.transpose() * d); + + if (δ > -a * σ) { + return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2)); + } else if (δ > -2 * a * σ) { + return 0.5 * k * pow(δ + 2 * a * σ, 2); + } else { + return 0; + } + }; +} + template <int D, class S> std::function<double(Spin<double, D, S>)> bCenter(double H) { return [H](Spin<double, D, S> s) -> double { return H * s.x.norm(); }; } -- cgit v1.2.3-70-g09d2