summaryrefslogtreecommitdiff
path: root/src/new_model.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/new_model.cpp')
-rw-r--r--src/new_model.cpp84
1 files changed, 48 insertions, 36 deletions
diff --git a/src/new_model.cpp b/src/new_model.cpp
index b96421d..326d530 100644
--- a/src/new_model.cpp
+++ b/src/new_model.cpp
@@ -56,19 +56,21 @@ double vsin(const vector_t<2, double>& v1, const vector_t<2, double>& v2) {
return v1[0] * v2[1] - v1[1] * v2[0];
}
+// these need to be global so that both wolff hooks have access to them.
+// TODO: investigate using pointers
double cos_xy;
vector_t<2, double> sin_xy;
class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i> {
private:
double a;
+ double K1;
- unsigned n;
int A;
public:
double totalA2;
- i_measurement(const ising::wolff::system<ising_t, ising_t, graph_i>& S, double a) : a(a) {
+ i_measurement(const ising::wolff::system<ising_t, ising_t, graph_i>& S, double a, double K1) : a(a), K1(K1) {
A = S.nv;
totalA2 = 0;
@@ -81,7 +83,7 @@ class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i
A -= s_new - s_old;
}
- sin_xy[v.prop.sublattice] += a * (s_new - s_old) * vsin(*(v.prop.back), *(v.prop.front));
+ sin_xy[v.prop.sublattice] += K1 * a * (s_new - s_old) * vsin(*(v.prop.back), *(v.prop.front));
cos_xy -= dE;
}
@@ -93,13 +95,13 @@ class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i
class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, vector_t<2, double>, graph_x> {
private:
- double Jxy;
- double Jxyn;
+ double K1y;
+ double K1yn;
double a;
const ising_t* si0;
public:
- 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) {
+ x_measurement(const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>, graph_x>& S, double K1y, double K1yn, double a, const ising_t *si0) : K1y(K1y), K1yn(K1yn), a(a), si0(si0) {
}
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 {
@@ -108,9 +110,9 @@ class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, vec
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;
if (e.prop.distance == 0) {
- sin_xy += (Jxy + (*si0) * a * *(e.prop.spin)) * dsin;
+ sin_xy += K1y * (1 + (*si0) * a * *(e.prop.spin)) * dsin;
} else if (e.prop.distance == 1) {
- sin_xy += Jxyn * dsin;
+ sin_xy += K1yn * dsin;
}
}
};
@@ -119,18 +121,19 @@ int main (int argc, char *argv[]) {
// set defaults
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;
+
+ 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 α = 1.5; // ising-xy coupling
int opt;
// take command line arguments
- while ((opt = getopt(argc, argv, "N:L:T:I:X:a:n:")) != -1) {
+ while ((opt = getopt(argc, argv, "N:L:T:J:K:a:n:")) != -1) {
switch (opt) {
case 'N': // number of steps
N = (unsigned)atof(optarg);
@@ -138,20 +141,20 @@ int main (int argc, char *argv[]) {
case 'L': // linear size
L = atoi(optarg);
break;
- case 'T': // temperature
+ case 'T':
T = atof(optarg);
break;
- case 'I': // temperature
- Ji = atof(optarg);
+ case 'J':
+ J = atof(optarg);
break;
- case 'X': // temperature
- Jx = atof(optarg);
+ case 'K':
+ K1 = atof(optarg);
break;
- case 'a': // temperature
- a = atof(optarg);
+ case 'n':
+ K2 = atof(optarg);
break;
- case 'n': // temperature
- Jxn = atof(optarg);
+ case 'a':
+ α = atof(optarg);
break;
default:
exit(EXIT_FAILURE);
@@ -165,16 +168,16 @@ int main (int argc, char *argv[]) {
// 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 -Ji;
+ return -J;
} else {
- return Ji;
+ return J;
}
};
// 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));
+ [K1, α] (const graph_i::vertex& v, const ising_t& s) -> double {
+ return K1 * α * s * (*(v.prop.back) * *(v.prop.front));
};
// initialize Ising diagonal-lattice graph
@@ -213,11 +216,11 @@ int main (int argc, char *argv[]) {
// 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 {
+ Zxy = [K1, K2, α, &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);
+ return K1 * (1 + si.s0 * α * *(e.prop.spin)) * (s1 * s2);
} else if (e.prop.distance == 1) {
- return Jxn * (s1 * s2);
+ return K1 * K2 * (s1 * s2);
} else {
return 0;
}
@@ -326,11 +329,11 @@ int main (int argc, char *argv[]) {
si.s[pow(L, D) + i].x = true;
}
- cos_xy = Jx * 2 * pow(L, 2) + Jxn * 2 * pow(L, 2);
+ cos_xy = K1 * 2 * pow(L, 2) + K2 * 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, α, K1);
+ x_measurement mxy(sxy, K1, K2, α, &(si.s0));
double total_cos = 0;
double total_sin2 = 0;
@@ -344,9 +347,18 @@ int main (int argc, char *argv[]) {
// std::cout << cos_xy / sxy.ne << " " << sin_xy / pow(sxy.ne, 2) << " " << E << "\n";
}
+ std::ifstream checkfile("out.dat");
+ bool already_exists = checkfile.good();
+ if (already_exists) {
+ checkfile.close();
+ }
+
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";
+ if (!already_exists) {
+ outfile << "N L J K1 K2 \\[Alpha] M \\[Rho] \\[Sigma]\n";
+ }
+ outfile << N << " " << L << " " << 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.close();
}