summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--dimers.cpp2
-rw-r--r--dimers_torus.cpp2
-rw-r--r--spheres.cpp2
-rw-r--r--spheres.hpp39
-rw-r--r--vector.hpp2
5 files changed, 24 insertions, 23 deletions
diff --git a/dimers.cpp b/dimers.cpp
index 4416333..d5cfdf4 100644
--- a/dimers.cpp
+++ b/dimers.cpp
@@ -85,7 +85,7 @@ int main(int argc, char* argv[]) {
file.open(filename);
Animation<double, D, Euclidean<double, D>, Dimer<double, D>> A(L, 750, argc, argv, wait);
- model sphere(1.0, Z, bCenter<D, Dimer<double, D>>(H));
+ model sphere(1, Z, bCenter<D, Dimer<double, D>>(H));
Rng rng;
diff --git a/dimers_torus.cpp b/dimers_torus.cpp
index 74cf47e..ec7c73a 100644
--- a/dimers_torus.cpp
+++ b/dimers_torus.cpp
@@ -86,7 +86,7 @@ int main(int argc, char* argv[]) {
sphere.s.resize(n);
- unsigned nx = floor(sqrt(n));
+ unsigned nx = ceil(sqrt(n));
for (unsigned i = 0; i < sphere.s.size(); i++) {
Spin<double, 2, Dimer<double, D>>* ss = new Spin<double, 2, Dimer<double, D>>();
ss->x = {(i / nx) * L / nx, (i % nx) * L / nx};
diff --git a/spheres.cpp b/spheres.cpp
index 9ca2036..4b2237d 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -75,7 +75,7 @@ int main(int argc, char* argv[]) {
sphere.s.resize(n);
- unsigned nx = floor(sqrt(n));
+ unsigned nx = ceil(sqrt(n));
for (unsigned i = 0; i < sphere.s.size(); i++) {
Spin<double, 2, double>* ss = new Spin<double, 2, double>();
ss->x = {(i / nx) * L / nx, (i % nx) * L / nx};
diff --git a/spheres.hpp b/spheres.hpp
index 3dc4a5e..60154f6 100644
--- a/spheres.hpp
+++ b/spheres.hpp
@@ -2,13 +2,8 @@
#include "space_wolff.hpp"
template <int D>
-std::function<double(const Spin<double, D, Radius>&, const Spin<double, D, Radius>&)>
-zSpheres(double a, double k) {
- return [a, k](const Spin<double, D, Radius>& s1, const Spin<double, D, Radius>& s2) -> double {
- Vector<double, D> d = s1.x - s2.x;
-
- double σ = s1.s + s2.s;
- double δ = σ - sqrt(d.transpose() * d);
+double softPotential(double a, double k, const Vector<double, D>& diff, double σ) {
+ double δ = σ - sqrt(diff.transpose() * diff);
if (δ > -a * σ) {
return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
@@ -17,6 +12,23 @@ zSpheres(double a, double k) {
} else {
return 0;
}
+}
+
+template <int D>
+double hardPotential(const Vector<double, D>& diff, double σ) {
+ if (pow(σ, 2) < diff.transpose() * diff) {
+ return 0;
+ } else {
+ return -std::numeric_limits<double>::infinity();
+ }
+}
+
+template <int D>
+std::function<double(const Spin<double, D, Radius>&, const Spin<double, D, Radius>&)>
+zSpheres(double a, double k) {
+ return [a, k](const Spin<double, D, Radius>& s1, const Spin<double, D, Radius>& s2) -> double {
+ Vector<double, D> d = s1.x - s2.x;
+ return softPotential(a, k, d, s1.s + s2.s);
};
}
@@ -24,18 +36,7 @@ template <int D>
std::function<double(const Spin<double, D, Radius>&, const Spin<double, D, Radius>&)>
zSpheresTorus(double L, double a, double k) {
return [L, a, k](const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
- Vector<double, D> d = diff(L, s1.x, s2.x);
-
- double σ = s1.s + s2.s;
- double δ = σ - sqrt(d.transpose() * d);
-
- if (δ > -a * σ) {
- return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
- } else if (δ > -2 * a * σ) {
- return 0.5 * k * pow(δ + 2 * a * σ, 2);
- } else {
- return 0;
- }
+ return softPotential(a, k, diff(L, s1.x, s2.x), s1.s + s2.s);
};
}
diff --git a/vector.hpp b/vector.hpp
index f642816..04a1a97 100644
--- a/vector.hpp
+++ b/vector.hpp
@@ -10,7 +10,7 @@ template <class T, int D> Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v
Vector<T, D> v;
for (unsigned i = 0; i < D; i++) {
- v(i) = std::abs(mod(v1(i), L) - mod(v2(i), L));
+ v(i) = mod(std::abs(v1(i) - v2(i)), L);
if (v(i) > L / 2) {
v(i) = L - v(i);
}