summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2019-01-29 12:30:17 -0500
committerJaron Kent-Dobias <jaron@kent-dobias.com>2019-01-29 12:30:17 -0500
commit0ac4b4869ef9a2d468afdcc2ceab296afe036279 (patch)
treea72c3ed51c2e1ead6623232275d2858fc6e98c85 /src
parent87312a12b896ff6e5b400a9877e5c2a2803857b2 (diff)
downloadcode-0ac4b4869ef9a2d468afdcc2ceab296afe036279.tar.gz
code-0ac4b4869ef9a2d468afdcc2ceab296afe036279.tar.bz2
code-0ac4b4869ef9a2d468afdcc2ceab296afe036279.zip
nnn using new wolff edge properties
Diffstat (limited to 'src')
-rw-r--r--src/new_model.cpp505
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();
}