summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--ciamarra.cpp75
-rw-r--r--distinguishable.cpp20
-rw-r--r--glass.hpp73
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);
diff --git a/glass.hpp b/glass.hpp
index 8ad0b52..1b1beb7 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -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;
- }
};