summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--distinguishable.cpp246
-rw-r--r--glass.hpp230
2 files changed, 240 insertions, 236 deletions
diff --git a/distinguishable.cpp b/distinguishable.cpp
index 4b66989..611bf50 100644
--- a/distinguishable.cpp
+++ b/distinguishable.cpp
@@ -1,5 +1,4 @@
#include <iostream>
-#include <queue>
#include "glass.hpp"
@@ -12,262 +11,39 @@ public:
DistinguishableState(unsigned t) { type = t; }
bool empty() const { return type == 0; }
+ bool operator==(const DistinguishableState& s) const { return type == s.type; }
void remove() { type = 0; }
};
-template <unsigned D> class System {
+template <unsigned D>
+class DistinguishableSystem : public FiniteEnergySystem<D, DistinguishableState> {
public:
- const unsigned L;
- unsigned N;
- std::vector<Vertex<D, DistinguishableState>> vertices;
std::vector<double> interaction;
- unsigned vectorToIndex(const Vector<D>& x) const {
- unsigned i = 0;
- for (unsigned d = 0; d < D; d++) {
- i += x[d] * iPow(L, d);
- }
- return i;
- }
-
- Vector<D> indexToVector(unsigned i) const {
- Vector<D> x;
- for (unsigned d = 0; d < D; d++) {
- x[d] = (i / iPow(L, d)) % L;
- }
- return x;
- }
-
- System(unsigned L, unsigned N, Rng& r) : L(L), N(N), vertices(iPow(L, D)), interaction(N * N) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
- vertices[i].position = indexToVector(i);
- vertices[i].initialPosition = vertices[i].position;
- vertices[i].adjacentEdges.reserve(2 * D);
- vertices[i].marked = false;
- }
-
- for (unsigned d = 0; d < D; d++) {
- Vector<D> Δx = Vector<D>::Zero();
- Δx[d] = 1;
- for (signed i = 0; i < iPow(L, D); i++) {
- unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1));
- vertices[i].adjacentEdges.push_back(HalfEdge<D, DistinguishableState>(vertices[j], Δx));
- vertices[j].adjacentEdges.push_back(HalfEdge<D, DistinguishableState>(vertices[i], -Δx));
- }
- }
-
+ DistinguishableSystem(unsigned L, unsigned N, Rng& r)
+ : FiniteEnergySystem<D, DistinguishableState>(L, N), interaction(N * N) {
for (double& V : interaction) {
V = r.uniform(-0.5, 0.5);
}
- for (unsigned i = 0; i < N; i++) {
- vertices[i].state.type = i + 1;
+ for (unsigned i = 0; i < this->N; i++) {
+ this->vertices[i].state.type = i + 1;
}
}
- double pairEnergy(const DistinguishableState& s1, const DistinguishableState& s2) const {
+ double pairEnergy(const DistinguishableState& s1, const DistinguishableState& s2) const override {
if (s1.empty() || s2.empty()) {
return 0;
} else {
if (s1.type < s2.type) {
- return interaction[(s1.type - 1) * N + (s2.type - 1)];
+ return interaction[(s1.type - 1) * this->N + (s2.type - 1)];
} else {
- return interaction[(s2.type - 1) * N + (s1.type - 1)];
+ return interaction[(s2.type - 1) * this->N + (s1.type - 1)];
}
}
}
-
- double siteEnergy(const Vertex<D, DistinguishableState>& v) const {
- double E = 0;
-
- for (const HalfEdge<D, DistinguishableState>& e : v.adjacentEdges) {
- E += pairEnergy(v.state, e.neighbor.state);
- }
-
- return E;
- }
-
- double energy() const {
- double E = 0;
-
- for (const Vertex<D, DistinguishableState>& v : vertices) {
- E += siteEnergy(v);
- }
-
- return E / 2;
- }
-
- bool trySwap(Vertex<D, DistinguishableState>& v1, Vertex<D, DistinguishableState>& v2, double T,
- Rng& r) {
- double E0 = siteEnergy(v1) + siteEnergy(v2);
- std::swap(v1.state, v2.state);
- double E1 = siteEnergy(v1) + siteEnergy(v2);
- double ΔE = E1 - E0;
-
- if (exp(-ΔE / T) > r.uniform(0.0, 1.0)) {
- /* Accept the swap. */
- std::swap(v1.initialPosition, v2.initialPosition);
- return true;
- } else {
- /* Revert the swap. */
- std::swap(v1.state, v2.state);
- return false;
- }
- }
-
- bool tryRandomMove(double T, Rng& r) {
- Vertex<D, DistinguishableState>& v1 = r.pick(vertices);
-
- if (v1.empty()) {
- return false;
- }
-
- Vertex<D, DistinguishableState>& v2 = (r.pick(v1.adjacentEdges)).neighbor;
-
- if (!v2.empty()) {
- return false;
- }
-
- return trySwap(v1, v2, T, r);
- }
-
- bool tryRandomSwap(double T, Rng& r) {
- Vertex<D, DistinguishableState>& v1 = r.pick(vertices);
- Vertex<D, DistinguishableState>& v2 = r.pick(vertices);
-
- if (v1.state.type != v2.state.type) {
- return trySwap(v1, v2, T, r);
- } else {
- return false;
- }
- }
-
- double density() const { return N / pow(L, D); }
-
- void sweepLocal(Rng& r) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
- tryRandomMove(r);
- }
- }
-
- void sweepSwap(Rng& r) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
- tryRandomSwap(r);
- }
- }
-
- unsigned flipCluster(const Transformation<D>& R, Vertex<D, DistinguishableState>& v0, double T,
- Rng& r, bool dry = false) {
- std::queue<std::array<std::reference_wrapper<Vertex<D, DistinguishableState>>, 2>> q;
- Vector<D> x0New = R.apply(v0.position);
- Vertex<D, DistinguishableState>& v0New = vertices[vectorToIndex(x0New)];
- q.push({v0, v0New});
-
- unsigned n = 0;
-
- while (!q.empty()) {
- auto [vR, vNewR] = q.front();
- q.pop();
-
- Vertex<D, DistinguishableState>& v = vR;
- Vertex<D, DistinguishableState>& vNew = vNewR;
-
- if (!v.marked && !vNew.marked) {
- v.marked = true;
- vNew.marked = true;
-
- for (HalfEdge<D, DistinguishableState>& e : v.adjacentEdges) {
- Vertex<D, DistinguishableState>& vn = e.neighbor;
-
- Vector<D> xnNew = R.apply(vn.position);
- Vertex<D, DistinguishableState>& vnNew = vertices[vectorToIndex(xnNew)];
-
- if (!vn.marked && !vnNew.marked) {
- double E0 = pairEnergy(v.state, vn.state) + pairEnergy(vNew.state, vnNew.state);
- double E1 = pairEnergy(vNew.state, vn.state) + pairEnergy(v.state, vnNew.state);
- double ΔE = E1 - E0;
-
- if (exp(-ΔE / T) < r.uniform(0.0, 1.0)) {
- q.push({vn, vnNew});
- }
- }
- }
-
- if (!dry) {
- std::swap(v.state, vNew.state);
- std::swap(v.initialPosition, vNew.initialPosition);
- }
-
- n += 1;
- }
- }
-
- return n;
- }
-
- unsigned wolff(const Transformation<D>& R, double T, Rng& r) {
- unsigned n = flipCluster(R, r.pick(vertices), T, r);
- for (Vertex<D, DistinguishableState>& v : vertices) {
- v.marked = false;
- }
- return n;
- }
-
- void swendsenWang(const Transformation<D>& R, double T, Rng& r) {
- for (Vertex<D, DistinguishableState>& v : vertices) {
- if (!v.marked) {
- bool dry = 0.5 < r.uniform(0.0, 1.0);
- flipCluster(R, v, T, r, dry);
- }
- }
-
- for (Vertex<D, DistinguishableState>& v : vertices) {
- v.marked = false;
- }
- }
-
- int overlap(const System<D>& s) const {
- int o = 0;
-
- for (unsigned i = 0; i < vertices.size(); i++) {
- // DistinguishableState<D> s2 =
- // orientation.apply(s.vertices[vectorToIndex(orientation.inverse().apply(indexToVector(i)))].state);
- // o += vertices[i].state.dot(s2);
- }
-
- return o;
- }
-
- void setInitialPosition() {
- for (Vertex<D, DistinguishableState>& v : vertices) {
- v.initialPosition = v.position;
- }
- }
-
- double selfIntermediateScattering() const {
- double F = 0;
-
- Vector<D> k1 = {M_PI, 0};
- Vector<D> k2 = {0, M_PI};
- for (const Vertex<D, DistinguishableState>& v : vertices) {
- if (!v.empty()) {
- F += cos(k1.dot(v.position - v.initialPosition));
- F += cos(k2.dot(v.position - v.initialPosition));
- }
- }
-
- return F / (2 * N);
- }
};
-template <unsigned D> Vector<D> randomVector(unsigned L, Rng& r) {
- Vector<D> x;
- for (unsigned i = 0; i < D; i++) {
- x[i] = r.uniform((unsigned)0, L - 1);
- }
- return x;
-}
-
int main() {
const unsigned D = 2;
unsigned L = 40;
@@ -278,7 +54,7 @@ int main() {
Rng r;
- System<D> s(L, N, r);
+ DistinguishableSystem<D> s(L, N, r);
std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
diff --git a/glass.hpp b/glass.hpp
index fbff7d6..64fe8b7 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -1,4 +1,5 @@
#include <list>
+#include <queue>
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/src/Core/CwiseNullaryOp.h>
@@ -11,6 +12,14 @@ template <int D> using Matrix = Eigen::Matrix<int, D, D>;
using Rng = randutils::random_generator<pcg32>;
+template <unsigned D> Vector<D> randomVector(unsigned L, Rng& r) {
+ Vector<D> x;
+ for (unsigned i = 0; i < D; i++) {
+ x[i] = r.uniform((unsigned)0, L - 1);
+ }
+ return x;
+}
+
int iPow(int x, unsigned p) {
if (p == 0)
return 1;
@@ -139,9 +148,10 @@ public:
}
};
-template <typename T> concept State = requires(T v) {
+template <typename T> concept State = requires(T v, const T& v2) {
{ v.empty() } -> std::same_as<bool>;
{v.remove()};
+ { v.operator==(v2) } -> std::same_as<bool>;
};
template <unsigned D, State S> class HalfEdge;
@@ -165,3 +175,221 @@ public:
HalfEdge(Vertex<D, S>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
};
+template <unsigned D, State S> class System {
+public:
+ const unsigned L;
+ unsigned N;
+ std::vector<Vertex<D, S>> vertices;
+
+ unsigned vectorToIndex(const Vector<D>& x) const {
+ unsigned i = 0;
+ for (unsigned d = 0; d < D; d++) {
+ i += x[d] * iPow(L, d);
+ }
+ return i;
+ }
+
+ Vector<D> indexToVector(unsigned i) const {
+ Vector<D> x;
+ for (unsigned d = 0; d < D; d++) {
+ x[d] = (i / iPow(L, d)) % L;
+ }
+ return x;
+ }
+
+ System(unsigned L, unsigned N) : L(L), N(N), vertices(iPow(L, D)) {
+ for (unsigned i = 0; i < iPow(L, D); i++) {
+ vertices[i].position = indexToVector(i);
+ vertices[i].initialPosition = vertices[i].position;
+ vertices[i].adjacentEdges.reserve(2 * D);
+ vertices[i].marked = false;
+ }
+
+ for (unsigned d = 0; d < D; d++) {
+ Vector<D> Δx = Vector<D>::Zero();
+ Δx[d] = 1;
+ for (signed i = 0; i < iPow(L, D); i++) {
+ unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1));
+ vertices[i].adjacentEdges.push_back(HalfEdge<D, S>(vertices[j], Δx));
+ vertices[j].adjacentEdges.push_back(HalfEdge<D, S>(vertices[i], -Δx));
+ }
+ }
+ }
+
+ double density() const { return N / pow(L, D); }
+
+ void setInitialPosition() {
+ for (Vertex<D, S>& v : vertices) {
+ v.initialPosition = v.position;
+ }
+ }
+
+ double selfIntermediateScattering() const {
+ double F = 0;
+
+ Vector<D> k1 = {M_PI, 0};
+ Vector<D> k2 = {0, M_PI};
+ for (const Vertex<D, S>& v : vertices) {
+ if (!v.empty()) {
+ F += cos(k1.dot(v.position - v.initialPosition));
+ F += cos(k2.dot(v.position - v.initialPosition));
+ }
+ }
+
+ return F / (2 * this->N);
+ }
+};
+
+template <unsigned D, State S> class FiniteEnergySystem : public System<D, S> {
+public:
+ using System<D, S>::System;
+
+ virtual double pairEnergy(const S&, const S&) const { return 0; }
+
+ double siteEnergy(const Vertex<D, S>& v) const {
+ double E = 0;
+
+ for (const HalfEdge<D, S>& e : v.adjacentEdges) {
+ E += pairEnergy(v.state, e.neighbor.state);
+ }
+
+ return E;
+ }
+
+ double energy() const {
+ double E = 0;
+
+ for (const Vertex<D, S>& v : this->vertices) {
+ E += siteEnergy(v);
+ }
+
+ return E / 2;
+ }
+
+ bool trySwap(Vertex<D, S>& v1, Vertex<D, S>& v2, double T, Rng& r) {
+ double E0 = siteEnergy(v1) + siteEnergy(v2);
+ std::swap(v1.state, v2.state);
+ double E1 = siteEnergy(v1) + siteEnergy(v2);
+ double ΔE = E1 - E0;
+
+ if (exp(-ΔE / T) > r.uniform(0.0, 1.0)) {
+ /* Accept the swap. */
+ std::swap(v1.initialPosition, v2.initialPosition);
+ return true;
+ } else {
+ /* Revert the swap. */
+ std::swap(v1.state, v2.state);
+ return false;
+ }
+ }
+
+ bool tryRandomMove(double T, Rng& r) {
+ Vertex<D, S>& v1 = r.pick(this->vertices);
+
+ if (v1.empty()) {
+ return false;
+ }
+
+ Vertex<D, S>& v2 = (r.pick(v1.adjacentEdges)).neighbor;
+
+ if (!v2.empty()) {
+ return false;
+ }
+
+ return trySwap(v1, v2, T, r);
+ }
+
+ bool tryRandomSwap(double T, Rng& r) {
+ Vertex<D, S>& v1 = r.pick(this->vertices);
+ Vertex<D, S>& v2 = r.pick(this->vertices);
+
+ if (v1.state != v2.state) {
+ return trySwap(v1, v2, T, r);
+ } else {
+ return false;
+ }
+ }
+
+ void sweepLocal(Rng& r) {
+ for (unsigned i = 0; i < iPow(this->L, D); i++) {
+ tryRandomMove(r);
+ }
+ }
+
+ void sweepSwap(Rng& r) {
+ for (unsigned i = 0; i < iPow(this->L, D); i++) {
+ tryRandomSwap(r);
+ }
+ }
+
+ unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, double T, Rng& r,
+ bool dry = false) {
+ std::queue<std::array<std::reference_wrapper<Vertex<D, S>>, 2>> q;
+ Vector<D> x0New = R.apply(v0.position);
+ Vertex<D, S>& v0New = this->vertices[this->vectorToIndex(x0New)];
+ q.push({v0, v0New});
+
+ unsigned n = 0;
+
+ while (!q.empty()) {
+ auto [vR, vNewR] = q.front();
+ q.pop();
+
+ Vertex<D, S>& v = vR;
+ Vertex<D, S>& vNew = vNewR;
+
+ if (!v.marked && !vNew.marked) {
+ v.marked = true;
+ vNew.marked = true;
+
+ for (HalfEdge<D, S>& e : v.adjacentEdges) {
+ Vertex<D, S>& vn = e.neighbor;
+
+ Vector<D> xnNew = R.apply(vn.position);
+ Vertex<D, S>& vnNew = this->vertices[this->vectorToIndex(xnNew)];
+
+ if (!vn.marked && !vnNew.marked) {
+ double E0 = pairEnergy(v.state, vn.state) + pairEnergy(vNew.state, vnNew.state);
+ double E1 = pairEnergy(vNew.state, vn.state) + pairEnergy(v.state, vnNew.state);
+ double ΔE = E1 - E0;
+
+ if (exp(-ΔE / T) < r.uniform(0.0, 1.0)) {
+ q.push({vn, vnNew});
+ }
+ }
+ }
+
+ if (!dry) {
+ std::swap(v.state, vNew.state);
+ std::swap(v.initialPosition, vNew.initialPosition);
+ }
+
+ n += 1;
+ }
+ }
+
+ return n;
+ }
+
+ unsigned wolff(const Transformation<D>& R, double T, Rng& r) {
+ unsigned n = flipCluster(R, r.pick(this->vertices), T, r);
+ for (Vertex<D, S>& v : this->vertices) {
+ v.marked = false;
+ }
+ return n;
+ }
+
+ void swendsenWang(const Transformation<D>& R, double T, Rng& r) {
+ for (Vertex<D, S>& v : this->vertices) {
+ if (!v.marked) {
+ bool dry = 0.5 < r.uniform(0.0, 1.0);
+ flipCluster(R, v, T, r, dry);
+ }
+ }
+
+ for (Vertex<D, S>& v : this->vertices) {
+ v.marked = false;
+ }
+ }
+};
+