summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/new_model.cpp375
1 files changed, 375 insertions, 0 deletions
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 <getopt.h>
+#include <iostream>
+
+#include "../include/randutils/randutils.hpp"
+
+#include <functional>
+#include <vector>
+#include <queue>
+
+#include <wolff/types.h>
+#include <wolff/graph.hpp>
+
+using namespace wolff;
+
+namespace ising {
+#define WOLFF_SITE_DEPENDENCE
+#include <wolff.hpp>
+#include <wolff/models/ising.hpp>
+#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 <wolff.hpp>
+#include <wolff/models/vector.hpp>
+#include <wolff/models/orthogonal.hpp>
+#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<ising_t, ising_t> {
+ 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<ising_t, ising_t>& 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<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 {
+ 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;
+ }
+
+ double avgM2() {
+ return totalM2 / ns;
+ }
+
+ double avgC() {
+ return totalC / nc;
+ }
+
+ double avgC2() {
+ return totalC2 / nc;
+ }
+};
+
+class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, 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<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;
+ }
+
+ 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 <double(const ising_t&, const ising_t&)> Zi = [=] (const ising_t& s1, const ising_t s2) -> double {
+ if (s1.x == s2.x) {
+ return -J;
+ } else {
+ return J;
+ }
+ };
+
+ 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);
+ }
+
+ 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<v_t> 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<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;
+
+ 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<orthogonal_t<2, double>, 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";
+}
+