From 30d0fee3be1b899e93c5af7cf9de585071bacd44 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 2 Jul 2021 17:10:44 +0200 Subject: Work on the Ciamarra model. --- ciamarra.cpp | 85 +++++++++++++++++++++++++------------ distinguishable.cpp | 16 +++++-- 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 class CiamarraSystem : public HardSystem> { + public: using HardSystem>::HardSystem; std::list>>> @@ -83,6 +84,38 @@ class CiamarraSystem : public HardSystem> { } } } + + /* 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(CiamarraSystem::L); + for (Vertex>& v : CiamarraSystem::vertices) { + v.initialPosition = 4 * v.position + v.state; + } + } + + double selfIntermediateScattering(const std::vector>& ms) const override { + double F = 0; + + std::array, D> ks; + for (unsigned i = 0; i < D; i++) { + ks[i].setZero(); + ks[i](i) = 1; + } + for (const Vertex>& v : CiamarraSystem::vertices) { + if (!v.empty()) { + Vector Δx = ((4 * CiamarraSystem::orientation.inverse().apply(v.position)) + CiamarraSystem::orientation.inverse().apply(v.state)) - v.initialPosition; + for (const Vector& 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 s0 = s; - std::vector> ms = generateTorusMatrices(); - 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 m = r.pick(ms); + Vertex>& v = r.pick(s.vertices); + unsigned n = s.wolff(Transformation(L, m, v.position - m * v.position), r); + if (i % 1 == 0) { + auto stop = std::chrono::high_resolution_clock::now(); + auto duration = duration_cast(stop - start); + std::cout << duration.count() << " " << s.selfIntermediateScattering(ms) << " " << n << " "; + } } - Matrix m = r.pick(ms); - Vertex>& v = r.pick(s.vertices); - unsigned nC = s.flipCluster(Transformation(L, m, v.position - m * v.position), v, r); - std::cerr << nC << std::endl; - for (Vertex>& 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(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 m = r.pick(ms); Vertex& v = r.pick(s.vertices); s.wolff(Transformation(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(stop - start); + std::cout << duration.count() << " " << s.selfIntermediateScattering(ms) << " "; + } } /* for (unsigned i = 0; i < 1e5; i++) { diff --git a/glass.hpp b/glass.hpp index fc7b186..8ad0b52 100644 --- a/glass.hpp +++ b/glass.hpp @@ -111,6 +111,8 @@ template std::vector> generateTorusMatrices() { return mats; } +template class CiamarraState; + template class Transformation { public: unsigned L; @@ -146,6 +148,11 @@ public: Transformation inverse() const { return Transformation(L, m.transpose(), -m.transpose() * v); } + + CiamarraState apply(const CiamarraState& s) const { + CiamarraState s2(m * s); + return s2; + } }; template concept State = requires(T v, const T& v2) { @@ -225,25 +232,31 @@ public: return iPow(L, D); } - void setInitialPosition() { + virtual void setInitialPosition() { + orientation = Transformation(L); for (Vertex& v : vertices) { v.initialPosition = v.position; } } - double selfIntermediateScattering() const { + virtual double selfIntermediateScattering(const std::vector>& ms) const { double F = 0; - Vector k1 = {M_PI, 0}; - Vector k2 = {0, M_PI}; + std::array, D> ks; + for (unsigned i = 0; i < D; i++) { + ks[i].setZero(); + ks[i](i) = 1; + } for (const Vertex& v : vertices) { if (!v.empty()) { - F += cos(k1.dot(v.position - v.initialPosition)); - F += cos(k2.dot(v.position - v.initialPosition)); + Vector Δx = v.position - orientation.apply(v.initialPosition); + for (const Vector& 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& 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& v : this->vertices) { v.marked = false; } @@ -431,6 +453,7 @@ public: bool tryRandomMove(Rng& r) { Vertex& v = r.pick(HardSystem::vertices); S oldState = v.state; + Vector oldPos = v.initialPosition; if (!tryDeletion(v)) { return false; @@ -440,6 +463,7 @@ public: for (HalfEdge& 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& v = r.pick(HardSystem::vertices); + S oldState = v.state; + Vector oldPos = v.initialPosition; + + if (!tryDeletion(v)) { + return false; + } + + Vertex& 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& v1, Vertex& 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 xNew = R.apply(v.position); Vertex& 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>> overlaps1 = overlaps(vNew, s, true); - std::list>> overlaps2 = overlaps(v, sNew, true); - overlaps1.splice(overlaps1.begin(), overlaps2); + if (s != vNew.state) { + v.marked = true; + vNew.marked = true; + + std::list>> overlaps1 = overlaps(vNew, s, true); + std::list>> overlaps2 = overlaps(v, sNew, true); + overlaps1.splice(overlaps1.begin(), overlaps2); - for (Vertex& vn : overlaps1) { - if (!vn.marked) { - q.push(vn); + for (Vertex& 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& 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& v : this->vertices) { + v.marked = false; + } + return n; + } + void swendsenWang(const Transformation& R, Rng& r) { for (Vertex& v : HardSystem::vertices) { if (!v.marked) { -- cgit v1.2.3-70-g09d2