From 87312a12b896ff6e5b400a9877e5c2a2803857b2 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sun, 21 Oct 2018 00:50:19 -0400 Subject: initialized git repo for work with Sam --- src/new_model.cpp | 375 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 375 insertions(+) create mode 100644 src/new_model.cpp (limited to 'src') diff --git a/src/new_model.cpp b/src/new_model.cpp new file mode 100644 index 0000000..96789d1 --- /dev/null +++ b/src/new_model.cpp @@ -0,0 +1,375 @@ + +#include +#include + +#include "../include/randutils/randutils.hpp" + +#include +#include +#include + +#include +#include + +using namespace wolff; + +namespace ising { +#define WOLFF_SITE_DEPENDENCE +#include +#include +#undef WOLFF_SITE_DEPEDENCE +} + +namespace xy { +#undef WOLFF_H +#undef WOLFF_SYSTEM_H +#undef WOLFF_CLUSTER_H +#undef WOLFF_MEASUREMENTS_H + +#define WOLFF_NO_FIELD +#define WOLFF_BOND_DEPENDENCE +#include +#include +#include +#undef WOLFF_NO_FIELD +#undef WOLFF_BOND_DEPENDENCE +} + +double E; + +using ising::wolff::ising_t; +using xy::wolff::vector_t; +using xy::wolff::orthogonal_t; + +class i_measurement : public ising::wolff::measurement { + private: + N_t nc; + N_t ns; + + int M; + int A; + v_t C; + + double totalaA; + double totalA2; + double totalaM; + double totalM2; + double totalC; + double totalC2; + + public: + i_measurement(const ising::wolff::system& S) { + nc = 0; + ns = 0; + M = 0; + A = S.nv; + + totalaA = 0; + totalA2 = 0; + totalaM = 0; + totalM2 = 0; + totalC = 0; + totalC2 = 0; + } + + void pre_cluster(N_t, N_t, const ising::wolff::system&, v_t, const ising_t&) override { + C = 0; + } + + void plain_bond_visited(const ising::wolff::system&, v_t, const ising_t&, v_t, double dE) override { + E += dE; + } + + void ghost_bond_visited(const ising::wolff::system& S, v_t i, const ising_t& s_old, const ising_t& s_new, double dE) override { + E += dE; + M += s_new - s_old; + + if (i >= S.nv / 2) { + A -= s_new - s_old; + } else { + A += s_new - s_old; + } + } + + void plain_site_transformed(const ising::wolff::system& S, v_t i, const ising_t& si_new) override { + C++; + } + + void post_cluster(N_t, N_t, const ising::wolff::system&) override { + totalaA += abs(A); + totalA2 += pow(A, 2); + totalaM += abs(M); + totalM2 += pow(M, 2); + totalC += C; + totalC2 += pow(C, 2); + + nc++; + + ns++; + } + + double avgaA() { + return totalaA / ns; + } + + double avgA2() { + return totalA2 / ns; + } + + double avgaM() { + return totalaM / ns; + } + + double avgM2() { + return totalM2 / ns; + } + + double avgC() { + return totalC / nc; + } + + double avgC2() { + return totalC2 / nc; + } +}; + +class x_measurement : public xy::wolff::measurement, vector_t<2, double>> { + private: + N_t ns; + N_t nc; + orthogonal_t<2, double> R; + + vector_t<2, double> M; + v_t C; + vector_t<2, double> CM; + + vector_t<2, double> totalaM; + vector_t<2, double> totalM2; + double totalM; + double totalCM; + double totalC; + double totalC2; + + public: + x_measurement(const xy::wolff::system, vector_t<2, double>>& S) { + ns = 0; + nc = 0; + M = S.nv * S.s[0]; + + totalaM.fill(0); + totalM2.fill(0); + totalC = 0; + totalC2 = 0; + totalM = 0; + totalCM = 0; + } + + void pre_cluster(N_t, N_t, const xy::wolff::system, vector_t<2, double>>&, v_t, const orthogonal_t<2, double>& r) override { + R = r; + C = 0; + CM.fill(0); + } + + void plain_bond_visited(const xy::wolff::system, vector_t<2, double>>&, v_t, const vector_t<2, double>&, v_t, double dE) override { + E += dE; + } + + void plain_site_transformed(const xy::wolff::system, vector_t<2, double>>& S, v_t i, const vector_t<2, double>& si_new) override { + C++; + CM += S.s[i]; + M += si_new - S.s[i]; + } + + void post_cluster(N_t, N_t, const xy::wolff::system, vector_t<2, double>>&) override { + totalaM[0] += fabs(M[0]); + totalaM[1] += fabs(M[1]); + totalM2[0] += pow(M[0], 2); + totalM2[1] += pow(M[1], 2); + totalM += (pow(M[0], 2) + pow(M[1], 2)); + totalCM += pow(CM[0] * R[0][0] + CM[1] * R[0][1], 2) / C; + totalC += C; + totalC2 += pow(C, 2); + + nc++; + + ns++; + } + + vector_t<2, double> avgaM() { + return totalaM / ns; + } + + vector_t<2, double> avgM2() { + return totalM2 / ns; + } + + double avgM() { + return totalM / ns; + } + + double avgCM() { + return totalCM / ns; + } + + double avgC() { + return totalC / nc; + } + + double avgC2() { + return totalC2 / nc; + } +}; + +using namespace wolff; + +int main (int argc, char *argv[]) { + + // set defaults + N_t N = (N_t)1e4; + const D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double J = 1.0; + double a = 1.5; + + int opt; + + // take command line arguments + while ((opt = getopt(argc, argv, "N:L:T:J:a:")) != -1) { + switch (opt) { + case 'N': // number of steps + N = (N_t)atof(optarg); + break; + case 'L': // linear size + L = atoi(optarg); + break; + case 'T': // temperature + T = atof(optarg); + break; + case 'J': // temperature + J = atof(optarg); + break; + case 'a': // temperature + a = atof(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + std::function Zi = [=] (const ising_t& s1, const ising_t s2) -> double { + if (s1.x == s2.x) { + return -J; + } else { + return J; + } + }; + + std::vector> *xy_spins; + + std::function Bi = [L, D, a, &xy_spins] (v_t v, const ising_t& s) -> double { + v_t xyv1 = v % (v_t)pow(L, D); + v_t xyv2; + if (v / (v_t)pow(L, D) == 0) { + xyv2 = L * (xyv1 / L) + ((xyv1 % L) + 1) % L; + } else { + xyv2 = L * (((xyv1 / L) + 1) % L) + (xyv1 % L); + } + + double value = a * ((xy_spins->at(xyv1)) * (xy_spins->at(xyv2))); + + if (s.x) { + return -value; + } else { + return value; + } + }; + + graph Gi; + + Gi.nv = 2 * pow(L, 2); + Gi.ne = 2 * Gi.nv; + + Gi.adjacency.resize(Gi.nv); + Gi.coordinate.resize(Gi.nv); + + for (std::vector adj_i : Gi.adjacency) { + adj_i.reserve(4); + } + + for (D_t i = 0; i < 2; i++) { + v_t sb = i * pow(L, 2); + + for (v_t j = 0; j < pow(L, 2); j++) { + v_t vc = sb + j; + + Gi.adjacency[vc].push_back(((i + 1) % 2) * pow(L, 2) + j); + Gi.adjacency[vc].push_back(((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * (j / (v_t)pow(L, 2 - 1)) + (j + 1 - 2 * (i % 2)) % L); + Gi.adjacency[vc].push_back(((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * ((L + (j/ (v_t)pow(L, 2 - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L); + Gi.adjacency[vc].push_back(((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * ((L + (j/ (v_t)pow(L, 2 - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L); + } + } + + ising::wolff::system si(Gi, T, Zi, Bi); + + std::function &, v_t, const vector_t<2, double>&)> Zxy = [L, D, a, &si] (v_t v1, const vector_t<2, double>& s1, v_t v2, const vector_t<2, double>& s2) -> double { + v_t iv; + if (v1 / L == v2 / L) { + v_t vs = v1 < v2 ? v1 : v2; + v_t vl = v1 < v2 ? v2 : v1; + + if (vl == L - 1 && vs == 0) { + iv = vl; + } else { + iv = vs; + } + } else { + v_t vs = v1 < v2 ? v1 : v2; + v_t vl = v1 < v2 ? v2 : v1; + + if (vl / L == L - 1 && vs / L == 0) { + iv = vl; + } else { + iv = vs; + } + } + + if (si.s[iv].x) { + return (1 + a) * (s1 * s2); + } else { + return (1 - a) * (s1 * s2); + } + }; + + graph Gxy(2, L); + + xy::wolff::system, vector_t<2, double>> sxy(Gxy, T, Zxy); + + xy_spins = &(sxy.s); + + for (v_t i = 0; i < pow(L, D); i++) { + si.s[pow(L, D) + i].x = true; + } + + E = - J * 4 * pow(L, D) - 2 * pow(L, D); + + i_measurement mi(si); + x_measurement mxy(sxy); + + double total_E; + double total_E2; + + for (N_t i = 0; i < N; i++) { + si.run_wolff(1, ising::wolff::gen_ising, mi, rng); + sxy.run_wolff(1, xy::wolff::generate_rotation_uniform<2>, mxy, rng); + total_E += E; + total_E2 += pow(E, 2); + } + + std::cout << mi.avgC() / si.nv << " " << mi.avgA2() / pow(si.nv, 2) << "\n"; + std::cout << mxy.avgM() / pow(sxy.nv, 1) << " " << 2 * mxy.avgCM() / pow(sxy.nv, 0) << "\n"; +} + -- cgit v1.2.3-70-g09d2