diff options
| -rw-r--r-- | ciamarra.cpp | 85 | ||||
| -rw-r--r-- | distinguishable.cpp | 16 | ||||
| -rw-r--r-- | glass.hpp | 119 |
3 files changed, 167 insertions, 53 deletions
diff --git a/ciamarra.cpp b/ciamarra.cpp index 7d7c8af..d1087a4 100644 --- a/ciamarra.cpp +++ b/ciamarra.cpp @@ -36,6 +36,7 @@ public: template <unsigned D> class CiamarraSystem : public HardSystem<D, CiamarraState<D>> { + public: using HardSystem<D, CiamarraState<D>>::HardSystem; std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>> @@ -83,6 +84,38 @@ class CiamarraSystem : public HardSystem<D, CiamarraState<D>> { } } } + + /* For the Ciamarra system, position within sites must be specially acconted + * for in the scattering function. We therefore expand the lattice and add + * the state when setting positions, so as later to be able to take particle + * sublattice positions into account. + * */ + void setInitialPosition() override { + CiamarraSystem::orientation = Transformation<D>(CiamarraSystem::L); + for (Vertex<D, CiamarraState<D>>& v : CiamarraSystem::vertices) { + v.initialPosition = 4 * v.position + v.state; + } + } + + double selfIntermediateScattering(const std::vector<Matrix<D>>& ms) const override { + double F = 0; + + std::array<Vector<D>, D> ks; + for (unsigned i = 0; i < D; i++) { + ks[i].setZero(); + ks[i](i) = 1; + } + for (const Vertex<D, CiamarraState<D>>& v : CiamarraSystem::vertices) { + if (!v.empty()) { + Vector<D> Δx = ((4 * CiamarraSystem::orientation.inverse().apply(v.position)) + CiamarraSystem::orientation.inverse().apply(v.state)) - v.initialPosition; + for (const Vector<D>& k : ks) { + F += cos((M_PI / 4) * k.dot(Δx)); + } + } + } + + return F / (D * CiamarraSystem::N); + } }; @@ -123,37 +156,35 @@ int main() { double z = exp(1 / 0.08); - while (s.density() < 0.818) { - s.sweepGrandCanonical(z, r); - } - - if (!s.compatible()) { - std::cerr << "Net compatible!" << std::endl; - return 1; - } - - std::cerr << "Found state with appropriate density." << std::endl; - - CiamarraSystem<D> s0 = s; - std::vector<Matrix<D>> ms = generateTorusMatrices<D>(); - unsigned n = 1; - for (unsigned i = 0; i < 1e5; i++) { - if (n < 20 * log(i + 1)) { - n++; - std::cout << i << " " << s.overlap(s0) << std::endl; + for (double ρ : {0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.76, 0.78, 0.8, 0.81, 0.815, 0.818, 0.821, 0.825}) { + while (s.density() < ρ) { + s.sweepGrandCanonical(z, r); + } + std::cout << s.density() << " "; + s.setInitialPosition(); + auto start = std::chrono::high_resolution_clock::now(); + for (unsigned i = 0; i < 1e4; i++) { + /* + for (unsigned j = 0; j < s.vertices.size(); j++) { + s.tryRandomNonlocalMove(r); + } + */ + Matrix<D> m = r.pick(ms); + Vertex<D, CiamarraState<D>>& v = r.pick(s.vertices); + unsigned n = s.wolff(Transformation<D>(L, m, v.position - m * v.position), r); + if (i % 1 == 0) { + auto stop = std::chrono::high_resolution_clock::now(); + auto duration = duration_cast<std::chrono::microseconds>(stop - start); + std::cout << duration.count() << " " << s.selfIntermediateScattering(ms) << " " << n << " "; + } } - Matrix<D> m = r.pick(ms); - Vertex<D, CiamarraState<D>>& v = r.pick(s.vertices); - unsigned nC = s.flipCluster(Transformation<D>(L, m, v.position - m * v.position), v, r); - std::cerr << nC << std::endl; - for (Vertex<D, CiamarraState<D>>& v : s.vertices) { - v.marked = false; + if (!s.compatible()) { + std::cerr << "Bad configuration" << std::endl; + return 1; } -// s.sweepLocal(r); -// s.sweepSwap(r); -// s.swendsenWang(Transformation<D>(L, ms, r), r); + std::cout << std::endl; } return 0; diff --git a/distinguishable.cpp b/distinguishable.cpp index 22f228a..b739927 100644 --- a/distinguishable.cpp +++ b/distinguishable.cpp @@ -83,13 +83,23 @@ int main() { */ /* */ - s.setInitialPosition(); std::cout << T << " "; - for (unsigned i = 0; i < 1e4; i++) { + s.setInitialPosition(); + auto start = std::chrono::high_resolution_clock::now(); + for (unsigned i = 0; i < 1e5; i++) { + for (unsigned j = 0; j < s.vertices.size(); j++) { + s.tryRandomSwap(T, r); + } + /* Matrix<D> m = r.pick(ms); Vertex<D, DistinguishableState>& v = r.pick(s.vertices); s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r); - std::cout << s.selfIntermediateScattering() << " "; + */ + if (i % 10 == 0) { + auto stop = std::chrono::high_resolution_clock::now(); + auto duration = duration_cast<std::chrono::microseconds>(stop - start); + std::cout << duration.count() << " " << s.selfIntermediateScattering(ms) << " "; + } } /* for (unsigned i = 0; i < 1e5; i++) { @@ -111,6 +111,8 @@ template <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() { return mats; } +template <unsigned D> class CiamarraState; + template <unsigned D> class Transformation { public: unsigned L; @@ -146,6 +148,11 @@ public: Transformation<D> inverse() const { return Transformation<D>(L, m.transpose(), -m.transpose() * v); } + + CiamarraState<D> apply(const CiamarraState<D>& s) const { + CiamarraState<D> s2(m * s); + return s2; + } }; template <typename T> concept State = requires(T v, const T& v2) { @@ -225,25 +232,31 @@ public: return iPow(L, D); } - void setInitialPosition() { + virtual void setInitialPosition() { + orientation = Transformation<D>(L); for (Vertex<D, S>& v : vertices) { v.initialPosition = v.position; } } - double selfIntermediateScattering() const { + virtual double selfIntermediateScattering(const std::vector<Matrix<D>>& ms) const { double F = 0; - Vector<D> k1 = {M_PI, 0}; - Vector<D> k2 = {0, M_PI}; + std::array<Vector<D>, D> ks; + for (unsigned i = 0; i < D; i++) { + ks[i].setZero(); + ks[i](i) = 1; + } 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)); + Vector<D> Δx = v.position - orientation.apply(v.initialPosition); + for (const Vector<D>& k : ks) { + F += cos(M_PI * k.dot(Δx)); + } } } - return F / (2 * this->N); + return F / (D * System::N); } }; @@ -371,7 +384,13 @@ public: std::swap(v.initialPosition, vNew.initialPosition); } - n += 1; + if (!v.empty()) { + n += 1; + } + + if (!vNew.empty()) { + n += 1; + } } } @@ -380,6 +399,9 @@ public: unsigned wolff(const Transformation<D>& R, double T, Rng& r) { unsigned n = flipCluster(R, r.pick(this->vertices), T, r); + if (n > SoftSystem::N / 2) { + SoftSystem::orientation = R.apply(SoftSystem::orientation); + } for (Vertex<D, S>& v : this->vertices) { v.marked = false; } @@ -431,6 +453,7 @@ public: bool tryRandomMove(Rng& r) { Vertex<D, S>& v = r.pick(HardSystem::vertices); S oldState = v.state; + Vector<D> oldPos = v.initialPosition; if (!tryDeletion(v)) { return false; @@ -440,6 +463,7 @@ public: for (HalfEdge<D, S>& e : v.adjacentEdges) { if (1 == e.Δx.dot(oldState)) { if (insert(e.neighbor, oldState.flip())) { + e.neighbor.initialPosition = oldPos; return true; } break; @@ -451,17 +475,46 @@ public: newState = S(r); } if (insert(v, newState)) { + v.initialPosition = oldPos; return true; } } + v.initialPosition = oldPos; + v.state = oldState; + HardSystem::N++; + return false; + } + + bool tryRandomNonlocalMove(Rng& r) { + Vertex<D, S>& v = r.pick(HardSystem::vertices); + S oldState = v.state; + Vector<D> oldPos = v.initialPosition; + + if (!tryDeletion(v)) { + return false; + } + + Vertex<D, S>& vNew = r.pick(HardSystem::vertices); + S newState(r); + + if (insert(vNew, newState)) { + vNew.initialPosition = oldPos; + return true; + } + + v.initialPosition = oldPos; v.state = oldState; HardSystem::N++; return false; } bool trySwap(Vertex<D, S>& v1, Vertex<D, S>& v2) { + if (v1.state == v2.state) { + return false; + } if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) { std::swap(v1.state, v2.state); + std::swap(v1.initialPosition, v2.initialPosition); return true; } else { return false; @@ -539,34 +592,54 @@ public: Vector<D> xNew = R.apply(v.position); Vertex<D, S>& vNew = HardSystem::vertices[HardSystem::vectorToIndex(xNew)]; - v.marked = true; - vNew.marked = true; - S s = R.apply(v.state); S sNew = R.apply(vNew.state); - std::list<std::reference_wrapper<Vertex<D, S>>> overlaps1 = overlaps(vNew, s, true); - std::list<std::reference_wrapper<Vertex<D, S>>> overlaps2 = overlaps(v, sNew, true); - overlaps1.splice(overlaps1.begin(), overlaps2); + if (s != vNew.state) { + v.marked = true; + vNew.marked = true; + + std::list<std::reference_wrapper<Vertex<D, S>>> overlaps1 = overlaps(vNew, s, true); + std::list<std::reference_wrapper<Vertex<D, S>>> overlaps2 = overlaps(v, sNew, true); + overlaps1.splice(overlaps1.begin(), overlaps2); - for (Vertex<D, S>& vn : overlaps1) { - if (!vn.marked) { - q.push(vn); + for (Vertex<D, S>& vn : overlaps1) { + if (!vn.marked) { + q.push(vn); + } } - } - if (!dry) { - v.state = sNew; - vNew.state = s; - } + if (!dry) { + v.state = sNew; + vNew.state = s; + std::swap(v.initialPosition, vNew.initialPosition); + } + + if (!v.empty()) { + n += 1; + } - n += 1; + if (!vNew.empty()) { + n += 1; + } + } } } return n; } + unsigned wolff(const Transformation<D>& R, Rng& r) { + unsigned n = flipCluster(R, r.pick(this->vertices), r); + if (n > HardSystem::N / 2) { + HardSystem::orientation = R.apply(HardSystem::orientation); + } + for (Vertex<D, S>& v : this->vertices) { + v.marked = false; + } + return n; + } + void swendsenWang(const Transformation<D>& R, Rng& r) { for (Vertex<D, S>& v : HardSystem::vertices) { if (!v.marked) { |