summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/new_model.cpp167
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();
}