diff options
| -rw-r--r-- | ciamarra.cpp | 75 | ||||
| -rw-r--r-- | distinguishable.cpp | 20 | ||||
| -rw-r--r-- | glass.hpp | 73 |
3 files changed, 69 insertions, 99 deletions
diff --git a/ciamarra.cpp b/ciamarra.cpp index d3acc38..1d4538f 100644 --- a/ciamarra.cpp +++ b/ciamarra.cpp @@ -84,38 +84,6 @@ 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); - } }; void print(const CiamarraSystem<2>& s) { @@ -161,27 +129,38 @@ int main() { while (s.density() < ρ) { s.sweepGrandCanonical(z, r); } - std::cout << s.density() << " "; - s.setInitialPosition(); + + std::cout << s.N << std::endl; + + CiamarraSystem<D> s0 = s; + double last_overlap = 1; auto start = std::chrono::high_resolution_clock::now(); - for (unsigned i = 0; i < 1e4; i++) { - /* + + while (last_overlap > 0.1) { for (unsigned j = 0; j < s.vertices.size(); j++) { - s.tryRandomNonlocalMove(r); + s.tryRandomMove(r); } - */ + + auto stop = std::chrono::high_resolution_clock::now(); + auto duration = duration_cast<std::chrono::microseconds>(stop - start); + last_overlap = ((double)s.overlap(s0) / s.N - s.density() / 6) / (1 - s.density() / 6); + std::cout << duration.count() << " " << last_overlap << " "; + } + std::cout << std::endl; + + CiamarraSystem<D> s1 = s; + last_overlap = 1; + auto start1 = std::chrono::high_resolution_clock::now(); + + while (last_overlap > 0.1) { 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 << " "; - } - } - if (!s.compatible()) { - std::cerr << "Bad configuration" << std::endl; - return 1; + s.wolff(Transformation<D>(L, m, v.position - m * v.position), r.pick(s.vertices)); + + auto stop = std::chrono::high_resolution_clock::now(); + auto duration = duration_cast<std::chrono::microseconds>(stop - start1); + last_overlap = ((double)s.overlap(s1) / s.N - s.density() / 6) / (1 - s.density() / 6); + std::cout << duration.count() << " " << last_overlap << " "; } std::cout << std::endl; } diff --git a/distinguishable.cpp b/distinguishable.cpp index 0aae0c1..de71c47 100644 --- a/distinguishable.cpp +++ b/distinguishable.cpp @@ -53,7 +53,7 @@ int main() { unsigned N = 1560; double Tmin = 0; double Tmax = 0.4; - double δT = 0.02; + double δT = 0.05; Rng r; @@ -86,23 +86,29 @@ int main() { /* */ std::cout << T << " "; + for (unsigned i = 0; i < 1e3; i++) { + for (unsigned j = 0; j < s.vertices.size(); j++) { + s.tryRandomSwap(T, r); + } + } s.setInitialPosition(); + DistinguishableSystem<D> s0 = s; auto start = std::chrono::high_resolution_clock::now(); - Quantity<double> energy(1e3); +// Quantity<double> energy(1e3); 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); - */ if (i % 10 == 0) { - energy.add(s.energy()); +// energy.add(s.energy()); auto stop = std::chrono::high_resolution_clock::now(); auto duration = duration_cast<std::chrono::microseconds>(stop - start); - std::cout /*<< duration.count() << " "*/ << s.selfIntermediateScattering(ms) << " "; + std::cout << duration.count() << " " << s.overlap(s0) << " "; } } /* @@ -114,12 +120,14 @@ int main() { } */ std::cout << std::endl; + /* std::vector<double> rho = energy.ρ(); for (const double& x : rho) { std::cout << x << " "; } std::cout << std::endl; std::cerr << T << " " << s.energy() / N << std::endl; + */ // s.sweepLocal(r); // s.sweepSwap(r); // s.swendsenWang(Transformation<D>(L, ms, r), r); @@ -153,6 +153,10 @@ public: CiamarraState<D> s2(m * s); return s2; } + + template <class A> A apply(const A& s) const { + return s; + }; }; template <typename T> concept State = requires(T v, const T& v2) { @@ -166,7 +170,6 @@ template <unsigned D, State S> class HalfEdge; template <unsigned D, State S> class Vertex { public: Vector<D> position; - Vector<D> initialPosition; S state; std::vector<HalfEdge<D, S>> adjacentEdges; bool marked; @@ -208,7 +211,6 @@ public: System(unsigned L) : L(L), N(0), vertices(iPow(L, D)), orientation(L) { 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; } @@ -232,31 +234,22 @@ public: return iPow(L, D); } - virtual void setInitialPosition() { - orientation = Transformation<D>(L); - for (Vertex<D, S>& v : vertices) { - v.initialPosition = v.position; - } - } + int overlap(const System<D, S>& s) const { + unsigned o = 0; + Transformation<D> oRel = orientation.apply(s.orientation.inverse()); - virtual double selfIntermediateScattering(const std::vector<Matrix<D>>& ms) const { - double F = 0; + for (unsigned i = 0; i < vertices.size(); i++) { + if (!vertices[i].empty()) { + S s1 = vertices[i].state; + S s2 = oRel.apply(s.vertices[vectorToIndex(oRel.inverse().apply(indexToVector(i)))].state); - 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()) { - Vector<D> Δx = v.position - orientation.apply(v.initialPosition); - for (const Vector<D>& k : ks) { - F += cos(M_PI * k.dot(Δx)); + if (s1 == s2) { + o += 1; } } } - return F / (D * System::N); + return o; } }; @@ -294,7 +287,6 @@ public: 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. */ @@ -381,7 +373,6 @@ public: if (!dry) { std::swap(v.state, vNew.state); - std::swap(v.initialPosition, vNew.initialPosition); } if (!v.empty()) { @@ -453,7 +444,6 @@ 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; @@ -463,7 +453,6 @@ 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; @@ -475,11 +464,9 @@ public: newState = S(r); } if (insert(v, newState)) { - v.initialPosition = oldPos; return true; } } - v.initialPosition = oldPos; v.state = oldState; HardSystem::N++; return false; @@ -488,7 +475,6 @@ public: 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; @@ -498,11 +484,9 @@ public: S newState(r); if (insert(vNew, newState)) { - vNew.initialPosition = oldPos; return true; } - v.initialPosition = oldPos; v.state = oldState; HardSystem::N++; return false; @@ -514,7 +498,6 @@ public: } 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; @@ -578,7 +561,7 @@ public: } } - unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, Rng& r, bool dry = false) { + unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, bool dry = false) { std::queue<std::reference_wrapper<Vertex<D, S>>> q; q.push(v0); @@ -612,7 +595,6 @@ public: if (!dry) { v.state = sNew; vNew.state = s; - std::swap(v.initialPosition, vNew.initialPosition); } if (!v.empty()) { @@ -629,8 +611,19 @@ public: return n; } + unsigned wolff(const Transformation<D>& R, Vertex<D, S>& v0) { + unsigned n = flipCluster(R, v0); + if (n > HardSystem::N / 2) { + HardSystem::orientation = R.apply(HardSystem::orientation); + } + for (Vertex<D, S>& v : this->vertices) { + v.marked = false; + } + return n; + } + unsigned wolff(const Transformation<D>& R, Rng& r) { - unsigned n = flipCluster(R, r.pick(this->vertices), r); + unsigned n = flipCluster(R, r.pick(this->vertices)); if (n > HardSystem::N / 2) { HardSystem::orientation = R.apply(HardSystem::orientation); } @@ -644,7 +637,7 @@ public: for (Vertex<D, S>& v : HardSystem::vertices) { if (!v.marked) { bool dry = 0.5 < r.uniform(0.0, 1.0); - unsigned n = flipCluster(R, v, r, dry); + unsigned n = flipCluster(R, v, dry); if (n > pow(HardSystem::L, D) / 4 && !dry) { HardSystem::orientation = R.apply(HardSystem::orientation); } @@ -656,14 +649,4 @@ public: } } - int overlap(const System<D, S>& s) const { - int o = 0; - - for (unsigned i = 0; i < HardSystem::vertices.size(); i++) { - S s2 = HardSystem::orientation.apply(s.vertices[HardSystem::vectorToIndex(HardSystem::orientation.inverse().apply(HardSystem::indexToVector(i)))].state); - o += HardSystem::vertices[i].state.dot(s2); - } - - return o; - } }; |