diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/new_model.cpp | 505 |
1 files changed, 241 insertions, 264 deletions
diff --git a/src/new_model.cpp b/src/new_model.cpp index 96789d1..b96421d 100644 --- a/src/new_model.cpp +++ b/src/new_model.cpp @@ -1,375 +1,352 @@ #include <getopt.h> #include <iostream> +#include <fstream> #include "../include/randutils/randutils.hpp" #include <functional> #include <vector> #include <queue> - -#include <wolff/types.h> -#include <wolff/graph.hpp> - -using namespace wolff; +#include <list> namespace ising { #define WOLFF_SITE_DEPENDENCE -#include <wolff.hpp> -#include <wolff/models/ising.hpp> +#include <wolff_models/ising.hpp> // ising includes wolff #undef WOLFF_SITE_DEPEDENCE + using namespace wolff; } 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 <wolff.hpp> -#include <wolff/models/vector.hpp> -#include <wolff/models/orthogonal.hpp> +#include <wolff_models/vector.hpp> +#include <wolff_models/orthogonal.hpp> // orthogonal includes wolff #undef WOLFF_NO_FIELD #undef WOLFF_BOND_DEPENDENCE + using namespace wolff; } -double E; +using ising::ising_t; +using xy::vector_t; +using xy::orthogonal_t; + +// at each ising site, we need to know which sublattice we're in and the xy +// spins to our left and right (top and bottom) +typedef struct _ising_site_props { + unsigned sublattice; + vector_t<2, double>* back; + vector_t<2, double>* front; +} ising_site_props; + +// at each xy bond, we need to know whether we're nearest or next-nearest, the +// displacement the bond makes, and the ising spin that lies on us (if any) +typedef struct _xy_bond_props { + unsigned distance; + vector_t<2, int> dNeighborSelf; + ising_t* spin; +} xy_bond_props; + +typedef ising::graph<ising_site_props> graph_i; +typedef xy::graph<std::tuple<>, xy_bond_props> graph_x; + +double vsin(const vector_t<2, double>& v1, const vector_t<2, double>& v2) { + return v1[0] * v2[1] - v1[1] * v2[0]; +} -using ising::wolff::ising_t; -using xy::wolff::vector_t; -using xy::wolff::orthogonal_t; +double cos_xy; +vector_t<2, double> sin_xy; -class i_measurement : public ising::wolff::measurement<ising_t, ising_t> { +class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i> { private: - N_t nc; - N_t ns; + double a; - int M; + unsigned n; int A; - v_t C; - - double totalaA; - double totalA2; - double totalaM; - double totalM2; - double totalC; - double totalC2; public: - i_measurement(const ising::wolff::system<ising_t, ising_t>& S) { - nc = 0; - ns = 0; - M = 0; + double totalA2; + i_measurement(const ising::wolff::system<ising_t, ising_t, graph_i>& S, double a) : a(a) { 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<ising_t, ising_t>&, v_t, const ising_t&) override { - C = 0; - } - - void plain_bond_visited(const ising::wolff::system<ising_t, ising_t>&, v_t, const ising_t&, v_t, double dE) override { - E += dE; - } - - void ghost_bond_visited(const ising::wolff::system<ising_t, ising_t>& 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 { + void ghost_bond_visited(const ising::wolff::system<ising_t, ising_t, graph_i>& S, const graph_i::vertex& v, const ising_t& s_old, const ising_t& s_new, double dE) override { + if (v.prop.sublattice == 0) { A += s_new - s_old; + } else { + A -= s_new - s_old; } - } - - void plain_site_transformed(const ising::wolff::system<ising_t, ising_t>& S, v_t i, const ising_t& si_new) override { - C++; - } - - void post_cluster(N_t, N_t, const ising::wolff::system<ising_t, ising_t>&) 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; - } + sin_xy[v.prop.sublattice] += a * (s_new - s_old) * vsin(*(v.prop.back), *(v.prop.front)); - double avgM2() { - return totalM2 / ns; + cos_xy -= dE; } - double avgC() { - return totalC / nc; - } - - double avgC2() { - return totalC2 / nc; + void post_cluster(unsigned, unsigned, const ising::wolff::system<ising_t, ising_t, graph_i>&) override { + totalA2 += pow(A, 2); } }; -class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, vector_t<2, double>> { +class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, vector_t<2, double>, graph_x> { 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; + double Jxy; + double Jxyn; + double a; + const ising_t* si0; public: - x_measurement(const xy::wolff::system<orthogonal_t<2, double>, 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<orthogonal_t<2, double>, 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<orthogonal_t<2, double>, 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<orthogonal_t<2, double>, 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<orthogonal_t<2, double>, 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; + 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) { } - double avgM() { - return totalM / ns; - } + 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 { + cos_xy -= dE; - double avgCM() { - return totalCM / ns; - } + 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; - double avgC() { - return totalC / nc; - } - - double avgC2() { - return totalC2 / nc; + if (e.prop.distance == 0) { + sin_xy += (Jxy + (*si0) * a * *(e.prop.spin)) * dsin; + } else if (e.prop.distance == 1) { + sin_xy += Jxyn * dsin; + } } }; -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; + 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; int opt; // take command line arguments - while ((opt = getopt(argc, argv, "N:L:T:J:a:")) != -1) { + while ((opt = getopt(argc, argv, "N:L:T:I:X:a:n:")) != -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); + case 'N': // number of steps + N = (unsigned)atof(optarg); + break; + case 'L': // linear size + L = atoi(optarg); + break; + case 'T': // temperature + T = atof(optarg); + break; + case 'I': // temperature + Ji = atof(optarg); + break; + case 'X': // temperature + Jx = atof(optarg); + break; + case 'a': // temperature + a = atof(optarg); + break; + case 'n': // temperature + Jxn = atof(optarg); + break; + default: + exit(EXIT_FAILURE); } } + // initialize random numbers randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; + // 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 -J; + return -Ji; } else { - return J; + return Ji; } }; - std::vector<vector_t<2, double>> *xy_spins; - - std::function <double(v_t, const ising_t&)> 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); - } + // 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)); + }; - double value = a * ((xy_spins->at(xyv1)) * (xy_spins->at(xyv2))); + // initialize Ising diagonal-lattice graph + graph_i Gi; - if (s.x) { - return -value; - } else { - return value; - } - }; + Gi.L = L; + Gi.D = 2; - graph Gi; - Gi.nv = 2 * pow(L, 2); Gi.ne = 2 * Gi.nv; - Gi.adjacency.resize(Gi.nv); - Gi.coordinate.resize(Gi.nv); + Gi.vertices.resize(Gi.nv); - for (std::vector<v_t> adj_i : Gi.adjacency) { - adj_i.reserve(4); - } + for (unsigned i = 0; i < 2; i++) { + unsigned sb = i * pow(L, 2); - for (D_t i = 0; i < 2; i++) { - v_t sb = i * pow(L, 2); + for (unsigned j = 0; j < pow(L, 2); j++) { + unsigned vc = sb + j; + Gi.vertices[vc].ind = vc; + Gi.vertices[vc].prop.sublattice = i; - for (v_t j = 0; j < pow(L, 2); j++) { - v_t vc = sb + j; + graph_i::halfedge e1(Gi.vertices[vc], Gi.vertices[((i + 1) % 2) * pow(L, 2) + j]); + graph_i::halfedge e2(Gi.vertices[vc], Gi.vertices[((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * (j / (unsigned)pow(L, 2 - 1)) + (j + 1 - 2 * (i % 2)) % L]); + graph_i::halfedge e3(Gi.vertices[vc], Gi.vertices[((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * ((L + (j/ (unsigned)pow(L, 2 - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L]); + graph_i::halfedge e4(Gi.vertices[vc], Gi.vertices[((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * ((L + (j/ (unsigned)pow(L, 2 - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L]); - 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); + Gi.vertices[vc].edges.push_back(e1); + Gi.vertices[vc].edges.push_back(e2); + Gi.vertices[vc].edges.push_back(e3); + Gi.vertices[vc].edges.push_back(e4); } } - - ising::wolff::system<ising_t, ising_t> si(Gi, T, Zi, Bi); - std::function <double(v_t, const vector_t<2, double>&, 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; + ising::wolff::system<ising_t, ising_t, graph_i> si(Gi, T, Zi, Bi); - if (vl == L - 1 && vs == 0) { - iv = vl; + // 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 { + if (e.prop.distance == 0) { + return (Jx + si.s0 * a * *(e.prop.spin)) * (s1 * s2); + } else if (e.prop.distance == 1) { + return Jxn * (s1 * s2); } else { - iv = vs; + return 0; } - } else { - v_t vs = v1 < v2 ? v1 : v2; - v_t vl = v1 < v2 ? v2 : v1; + }; + + // initialize square-lattice xy graph + graph_x Gxy(2, L); + + // assign edge properties to the XY nearest-neighbor bonds + for (graph_x::vertex& v : Gxy.vertices) { + for (graph_x::halfedge& e : v.edges) { + e.prop.distance = 0; + + int v1 = e.self.ind; + int v2 = e.neighbor.ind; + + e.prop.dNeighborSelf[0] = (v2 % L) - (v1 % L); + e.prop.dNeighborSelf[1] = v2 / L - v1 / L; - if (vl / L == L - 1 && vs / L == 0) { - iv = vl; + for (unsigned i = 0; i < 2; i++) { + if (e.prop.dNeighborSelf[i] == L - 1) { + e.prop.dNeighborSelf[i] = -1; + } + if (e.prop.dNeighborSelf[i] == -(L - 1)) { + e.prop.dNeighborSelf[i] = 1; + } + } + + unsigned vs = v1 < v2 ? v1 : v2; + unsigned vl = v1 < v2 ? v2 : v1; + + if (v1 / L == v2 / L) { + if (vl % L == L - 1 && vs % L == 0) { + e.prop.spin = &(si.s[vl]); + } else { + e.prop.spin = &(si.s[vs]); + } } else { - iv = vs; + if (vl / L == L - 1 && vs / L == 0) { + e.prop.spin = &(si.s[pow(L, 2) + vl]); + } else { + e.prop.spin = &(si.s[pow(L, 2) + vs]); + } } } + } - if (si.s[iv].x) { - return (1 + a) * (s1 * s2); - } else { - return (1 - a) * (s1 * s2); - } - }; + // add the next-nearest-neighbor bonds + for (graph_x::vertex &v : Gxy.vertices) { + int i = v.ind / L; + int j = v.ind % L; + + graph_x::halfedge e1(v, Gxy.vertices[((i - 1) % L) * L + (j + 1) % L]); + graph_x::halfedge e2(v, Gxy.vertices[((i - 1) % L) * L + (j - 1) % L]); + graph_x::halfedge e3(v, Gxy.vertices[((i + 1) % L) * L + (j + 1) % L]); + graph_x::halfedge e4(v, Gxy.vertices[((i + 1) % L) * L + (j - 1) % L]); + + e1.prop.distance = 1; + e2.prop.distance = 1; + e3.prop.distance = 1; + e4.prop.distance = 1; + + e1.prop.dNeighborSelf[0] = 1; + e1.prop.dNeighborSelf[1] = -1; + + e2.prop.dNeighborSelf[0] = -1; + e2.prop.dNeighborSelf[1] = -1; - graph Gxy(2, L); + e3.prop.dNeighborSelf[0] = 1; + e3.prop.dNeighborSelf[1] = 1; - xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>> sxy(Gxy, T, Zxy); + e4.prop.dNeighborSelf[0] = -1; + e4.prop.dNeighborSelf[1] = 1; - xy_spins = &(sxy.s); + v.edges.push_back(e1); + v.edges.push_back(e2); + v.edges.push_back(e3); + v.edges.push_back(e4); + } + + xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>, graph_x> sxy(Gxy, T, Zxy); + + // now that the xy spins have been initialized, go back to the ising graph + // and add a reference to the ones on each bond + for (unsigned i = 0; i < 2; i++) { + unsigned sb = i * pow(L, 2); - for (v_t i = 0; i < pow(L, D); i++) { + for (unsigned j = 0; j < pow(L, 2); j++) { + unsigned vc = sb + j; + + unsigned xyv1 = j; + unsigned xyv2; + if (i == 0) { + xyv2 = L * (xyv1 / L) + ((xyv1 % L) + 1) % L; + } else { + xyv2 = L * (((xyv1 / L) + 1) % L) + (xyv1 % L); + } + + si.G.vertices[vc].prop.back = &(sxy.s[xyv1]); + si.G.vertices[vc].prop.front = &(sxy.s[xyv2]); + } + } + + // start in the ground state of the Ising model + for (unsigned 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); + cos_xy = Jx * 2 * pow(L, 2) + Jxn * 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); - x_measurement mxy(sxy); + double total_cos = 0; + double total_sin2 = 0; - double total_E; - double total_E2; + for (unsigned i = 0; i < N; i++) { + si.run_wolff(1, ising::wolff::gen_ising<graph_i>, mi, rng); + sxy.run_wolff(1, xy::wolff::generate_rotation_uniform<2, graph_x>, mxy, rng); + total_cos += cos_xy; + total_sin2 += sin_xy * sin_xy; - 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 << cos_xy / sxy.ne << " " << sin_xy / pow(sxy.ne, 2) << " " << E << "\n"; } - 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"; + 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"; + outfile.close(); } |