diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/new_model.cpp | 167 |
1 files changed, 151 insertions, 16 deletions
diff --git a/src/new_model.cpp b/src/new_model.cpp index 4393a3d..4f473da 100644 --- a/src/new_model.cpp +++ b/src/new_model.cpp @@ -62,19 +62,129 @@ double vsin(const vector_t<2, double>& v1, const vector_t<2, double>& v2) { double cos_xy; vector_t<2, double> sin_xy; +std::vector<double> ρ(unsigned N, const std::vector<double>& C, double avg) { + double C0 = C.front() / N; + double avg2 = pow(avg / N, 2); + + std::vector<double> ρtmp; + + for (double Ct : C) { + ρtmp.push_back((Ct / N - avg2) / (C0 - avg2)); + } + + return ρtmp; +} + +double τ(unsigned n, const std::vector<double>& C, double avg) { + double τtmp = 0.5; + unsigned M = 1; + double c = 8.0; + + std::vector<double> ρ_tmp = ρ(n, C, avg); + + while (c * τtmp > M && M < C.size()) { + τtmp += ρ_tmp[M]; + M++; + } + + return τtmp; +} + +class quantity { + private: + double total; + std::vector<double> C; + + public: + unsigned n; + std::list<double> hist; + quantity(unsigned lag) : C(lag) { + n = 0; + total = 0; + } + + void add(double x) { + hist.push_front(x); + if (hist.size() > C.size()) { + hist.pop_back(); + unsigned t = 0; + for (double a : hist) { + C[t] += a * x; + t++; + } + total += x; + n++; + } + } + + double avg() { + return total / n; + } + + double σ() { + return 2.0 / n * τ(n, C, total) * (C[0] / n - pow(this->avg(), 2)); + } + + double serr() { + return sqrt(this->σ()); + } +}; + +class quantity_pair { + private: + std::vector<double> Cab; + double total; + double T; + quantity a; + quantity b; + + public: + quantity_pair(unsigned lag, double T) : a(lag), b(lag), Cab(lag), T(T) {}; + + void add(double x, double y) { + a.add(x); + b.add(y); + + if (a.hist.size() == Cab.size()) { + unsigned t = 0; + auto a_it = a.hist.begin(); + auto b_it = b.hist.begin(); + while (a_it != a.hist.end()) { + Cab[t] += (*a_it) * (*b_it); + t++; + a_it++; + b_it++; + } + } + } + + double cov() { + return 1.0 / a.n * τ(a.n, Cab, sqrt(a.avg() * b.avg()) * a.n) * (Cab[0] / a.n - a.avg() * b.avg()); + } + + double avg() { + return a.avg() - b.avg() / T; + } + + double serr() { + return sqrt(a.σ() + b.σ() / pow(T, 2) - 2 * this->cov() / T); + } +}; + class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i> { private: double a; double K1; - + unsigned w; int A; + unsigned N; public: - double totalA2; - i_measurement(const ising::wolff::system<ising_t, ising_t, graph_i>& S, double a, double K1) : a(a), K1(K1) { - A = S.nv; + quantity A2; - totalA2 = 0; + i_measurement(const ising::wolff::system<ising_t, ising_t, graph_i>& S, double a, double K1, unsigned nh, unsigned w) : a(a), w(w), K1(K1), A2(nh) { + A = S.nv; + N = 0; } 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 { @@ -90,7 +200,10 @@ class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i } void post_cluster(unsigned, unsigned, const ising::wolff::system<ising_t, ising_t, graph_i>&) override { - totalA2 += pow(A, 2); + if (N > w) { + A2.add(pow(A, 2)); + } + N++; } }; @@ -122,20 +235,24 @@ class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, vec int main (int argc, char *argv[]) { // set defaults - unsigned N = (unsigned)1e4; + unsigned N = (unsigned)1e4; // number of steps between precision checks + unsigned w = (unsigned)1e3; // steps to skip before data collection begins + unsigned nh = 1e2; // lag to keep in autocorrelation functions + double εmag = 1e-1; // precision goal in magnetization + double εstiff = 1e-1; // precision goal in spin stiffness const unsigned D = 2; // we're always in 2D unsigned L = 128; // system size double T = 1.0; // temperature double J = 1.0; // ising nn coupling double K1 = 1.0; // xy nn coupling - double K2 = 0.0; // xy nnn coupling + double K2 = 0.0; // xy nnn coupling double α = 1.5; // ising-xy coupling int opt; // take command line arguments - while ((opt = getopt(argc, argv, "N:L:T:J:K:a:n:")) != -1) { + while ((opt = getopt(argc, argv, "N:L:T:J:K:a:n:e:f:")) != -1) { switch (opt) { case 'N': // number of steps N = (unsigned)atof(optarg); @@ -158,6 +275,12 @@ int main (int argc, char *argv[]) { case 'a': α = atof(optarg); break; + case 'e': + εmag = atof(optarg); + break; + case 'f': + εstiff = atof(optarg); + break; default: exit(EXIT_FAILURE); } @@ -334,17 +457,29 @@ int main (int argc, char *argv[]) { cos_xy = K1 * 2 * pow(L, 2) + K1 * K2 * 4 * pow(L, 2); sin_xy.fill(0); - i_measurement mi(si, α, K1); + i_measurement mi(si, α, K1, nh, w); x_measurement mxy(sxy, K1, K2, α, &(si.s0)); - double total_cos = 0; - double total_sin2 = 0; + quantity_pair stiffness(nh, T); + + unsigned nn = 0; - for (unsigned i = 0; i < N; i++) { + while (true) { 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; + if (nn > w) { + stiffness.add(cos_xy, sin_xy * sin_xy); + + if (nn % N == 0) { + double err = mi.A2.serr(); + double val2 = stiffness.avg(); + double err2 = stiffness.serr(); + if (err / mi.A2.avg() <= εmag && err2 / val2 < εstiff) { + break; + } + } + } + nn++; } std::ifstream checkfile("out.dat"); @@ -358,7 +493,7 @@ int main (int argc, char *argv[]) { if (!already_exists) { outfile << "N L T J K1 K2 \\[Alpha] M \\[Rho] \\[Sigma]\n"; } - outfile << N << " " << L << " " << T << " " << J << " " << K1 << " " << K2 << " " << α << " " << 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 << N << " " << L << " " << T << " " << J << " " << K1 << " " << K2 << " " << α << " " << mi.A2.avg() / pow(si.nv, 2) << " " << mi.A2.serr() / pow(si.nv, 2) << " " << stiffness.avg() / sxy.ne << " " << stiffness.serr() / sxy.ne << "\n"; outfile.close(); } |