Add heatbath update for Heisenberg model & use it

frustmag/config_frustmag_sim.hpp

@@ -27,6 +27,7 @@
 #include "lattice/kagome.hpp"
 #include "lattice/dice.hpp"
 #include "update/single_flip.hpp"
+#include "update/heatbath.hpp"
 #include "update/overrelaxation.hpp"
 #include "update/mux.hpp"
 
@@ -96,7 +97,7 @@ template <template <typename> typename Lat>
 struct update_t<hamiltonian::heisenberg<Lat>> {
     template <typename LatticeH>
     using type = update::muxer<
-        update::single_flip
+        update::heatbath
         , update::overrelaxation
         , update::global_trafo
         , update::parallel_tempering

frustmag/hamiltonian/heisenberg.hpp

@@ -20,6 +20,7 @@
 #include "phase_space_point.hpp"
 #include "site/spin_O3.hpp"
 #include "update/single_flip.hpp"
+#include "update/heatbath.hpp"
 #include "update/overrelaxation.hpp"
 #include "update/global_trafo.hpp"
 #include "update/parallel_tempering.hpp"
@@ -96,6 +97,15 @@ struct heisenberg {
         return lattice_;
     }
 
+    auto local_field(site_iterator site_it) const {
+        auto nn = lattice().nearest_neighbors(site_it);
+        return std::accumulate(nn.begin(), nn.end(),
+           Eigen::Vector3d{Eigen::Vector3d::Zero()},
+           [] (auto const& a, auto const& b) {
+               return a + *b;
+           });
+    }
+
     template <typename RNG>
     // requires UniformRandomBitGenerator<RNG>
     bool metropolis(update::single_flip_proposal<lattice_type> const& p,
@@ -122,12 +132,7 @@ struct heisenberg {
     bool metropolis(update::overrelaxation_proposal<lattice_type> const& p,
                     RNG &)
     {
-        auto nn = lattice().nearest_neighbors(p.site_it);
-        auto n = std::accumulate(nn.begin(), nn.end(),
-                                 Eigen::Vector3d{Eigen::Vector3d::Zero()},
-                                 [] (auto const& a, auto const& b) {
-                                     return a + *b;
-                                 });
+        auto n = local_field(p.site_it);
         *(p.site_it) = site_state_type{2. * p.site_it->dot(n) / n.squaredNorm() * n
                                        - *(p.site_it)};
         return true;
@@ -164,6 +169,26 @@ struct heisenberg {
         return true;
     }
 
+    template <typename RNG>
+    // requires UniformRandomBitGenerator<RNG>
+    bool metropolis(update::heatbath_proposal<lattice_type> const& p,
+                    RNG & rng)
+    {
+        site_state_type m = local_field(p.site_it);
+        double alpha = beta * m.norm();
+        double phi_p = std::uniform_real_distribution<double>{0., 2. * M_PI}(rng);
+        double u = std::uniform_real_distribution<double>{}(rng);
+        double cos_theta_p = -1 - std::log(1 - u * (1 - std::exp(-2*alpha))) / alpha;
+        double sin_theta_p = std::sqrt(1 - std::pow(cos_theta_p, 2));
+        site_state_type s_p;
+        s_p << sin_theta_p * std::cos(phi_p),
+               sin_theta_p * std::sin(phi_p),
+               cos_theta_p;
+        *p.site_it = m.relative_spin(s_p);
+        return true;
+    }
+
+
     double log_weight(phase_point const& other) const {
         return (ppoint.temp - other.temp) * energy();
     }

frustmag/site/spin_O3.hpp

@@ -42,19 +42,19 @@ struct spin_O3 : Eigen::Vector3d {
     spin_O3() = default;
 
     spin_O3(Eigen::Vector3d const& other) : Eigen::Vector3d{other} {}
-    
-    template <typename RNG>
-    spin_O3 flipped(RNG & rng, double cos_theta_0 = -1) const {
-        double cos_theta = (*this)[2];
+
+    spin_O3 relative_spin(spin_O3 const& other) const {
+        double norm = this->norm();
+        double cos_theta = (*this)[2] / norm;
         double sin_theta = std::sqrt(1 - std::pow(cos_theta, 2));
-        double cos_phi = (*this)[0] / sin_theta;
-        double sin_phi = (*this)[1] / sin_theta;
+        double cos_phi = (*this)[0] / norm / sin_theta;
+        double sin_phi = (*this)[1] / norm / sin_theta;
         Eigen::Matrix3d R;
         R <<
             cos_theta * cos_phi, -sin_phi, sin_theta * cos_phi,
             cos_theta * sin_phi, cos_phi, sin_theta * sin_phi,
             -sin_theta, 0, cos_theta;
-        spin_O3 ret{R * random(rng, cos_theta_0)};
+        spin_O3 ret{R * other};
 
         // important: renormalize to avoid exponential growth of rounding errors
         ret /= ret.norm();
@@ -62,6 +62,11 @@ struct spin_O3 : Eigen::Vector3d {
         return ret;
     }
 
+    template <typename RNG>
+    spin_O3 flipped(RNG & rng, double cos_theta_0 = -1) const {
+        return relative_spin(random(rng, cos_theta_0));
+    }
+
     static const size_t size = 3;
 
     template <typename OutputIterator>

frustmag/update/heatbath.hpp

+// SVM Order Parameters for Hidden Spin Order
+// Copyright (C) 2018-2019  Jonas Greitemann, Ke Liu, and Lode Pollet
+
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#pragma once
+
+#include "concepts.hpp"
+
+#include <alps/params.hpp>
+
+namespace update {
+
+template <typename Lat>
+struct heatbath_proposal {
+#ifdef USE_CONCEPTS
+    static_assert(Lattice<Lat>, "Lat is not a Lattice");
+#endif
+    using site_iterator = typename Lat::iterator;
+
+    site_iterator site_it;
+};
+
+template <typename LatticeH>
+struct heatbath {
+#ifdef USE_CONCEPTS
+    static_assert(LatticeHamiltonian<LatticeH>,
+                  "LatticeH is not a LatticeHamiltonian");
+#endif
+private:
+    using lattice_type = typename LatticeH::lattice_type;
+    using site_iterator = typename lattice_type::iterator;
+    using site_state_type = typename lattice_type::value_type;
+public:
+    using hamiltonian_type = LatticeH;
+    using acceptance_type = std::array<double, 1>;
+
+    static void define_parameters(alps::params & parameters) {}
+
+    template <typename... Args>
+    heatbath(alps::params const&, Args &&...) {}
+
+    template <typename RNG>
+    // requires SiteState<site_state_type, RNG>
+    acceptance_type update(LatticeH & hamiltonian, RNG & rng) {
+        auto site_it = hamiltonian.lattice().begin();
+        for (; site_it != hamiltonian.lattice().end(); ++site_it) {
+            hamiltonian.metropolis(heatbath_proposal<lattice_type>{site_it},
+                                   rng);
+        }
+        return {1};
+    }
+};
+
+}