diff --git a/src/hamiltonian/Heisenberg.h b/src/hamiltonian/Heisenberg.h index cb46ca49631942f7f45a30a7e4d22ae3688b8804..3c6794deeb58f764fed9bd8a04ede6d9a157ff95 100644 --- a/src/hamiltonian/Heisenberg.h +++ b/src/hamiltonian/Heisenberg.h @@ -95,11 +95,11 @@ class Heisenberg // ---- Initializer ---- template - void initstatespace(StateSpace& statespace, Lattice& grid, RNG& rng) const + void initstatespace(StateSpace& statespace, Lattice& grid, RNG& rng) const noexcept { for(decltype(grid.size()) i = 0; i < grid.size(); ++i) { - statespace[i] = rnddir(rng); + set_to_rnddir(rng, statespace[i]); } } }; @@ -156,7 +156,7 @@ namespace MARQOV * @param rng reference to the random number generator */ template - void draw(RNG& rng) {rdir = rnddir(rng);} + void draw(RNG& rng) noexcept {set_to_rnddir(rng, rdir);} /** Computes the Wolff coupling when attempting to add a spin to the cluster diff --git a/src/hamiltonian/Phi4.h b/src/hamiltonian/Phi4.h index 4aab49bf00c0129d4131943b8777d249a9123a8f..3517e43ebd7d29511db9a71d6b43741171838eb1 100644 --- a/src/hamiltonian/Phi4.h +++ b/src/hamiltonian/Phi4.h @@ -30,8 +30,9 @@ class Phi4_Initializer_Radial { double amp = 0.5; // Amplitude (TODO: make me a class parameter) - auto newdir = rnddir(rng); - auto nsv = osv + mult(amp,newdir); + StateVector newdir; + set_to_rnddir(rng, newdir); + auto nsv = osv + mult(amp, newdir); return nsv; }; @@ -225,7 +226,7 @@ namespace MARQOV * @param rng reference to the random number generator */ template - void draw(RNG& rng, StateVector& sv) {rdir = rnddir(rng);} + void draw(RNG& rng, StateVector& sv) noexcept {set_to_rnddir(rng, rdir);} /** Computes the Wolff coupling when attempting to add a spin to the cluster diff --git a/src/hamiltonian/util/initializers.h b/src/hamiltonian/util/initializers.h index bca0f197750f65af90c7cb13a05a28188abad93c..c8ac3667b3fd85920219f00c19a7b420ade11ce8 100644 --- a/src/hamiltonian/util/initializers.h +++ b/src/hamiltonian/util/initializers.h @@ -55,24 +55,25 @@ public: } }; -template -class SVInitializer, typename std::enable_if< -std::is_floating_point::value>::type> +template +class SVInitializer::value>::type> { - typedef std::array StateVector; + using StateVector = T; public: - /** Flip the Ising spin. + /** Flip the Spin. * * @tparam RNGCache the Random Number Generator. * - * @param svold the old state vector. - * @param rng a reference to the RNG of the Monte Carlo simulation. - * @return the new state vector with flipped spin + * @param svold The old state vector. + * @param rng A reference to the RNG of the Monte Carlo simulation. + * @return The new state vector with flipped spin */ template static StateVector newsv(const StateVector& svold, RNGCache& rng) { - return rnddir(rng); + StateVector retval; + set_to_rnddir(rng, retval); + return retval; } }; @@ -86,11 +87,11 @@ class Spin1_Initializer public: /** Draw a new state vector based on the old one. * - * @tparam RNG the Random Number Generator. + * @tparam RNG The Random Number Generator. * - * @param svold the old state vector. - * @param rng a reference to the RNG of the Monte Carlo simulation. - * @return the new state vector + * @param svold The old state vector. + * @param rng A reference to the RNG of the Monte Carlo simulation. + * @return The new state vector */ template static StateVector newsv(const StateVector& svold, RNGCache& rng) diff --git a/src/hamiltonian/util/randomdir.h b/src/hamiltonian/util/randomdir.h index ff4c80e215b7de9df4e8de09cf24858e022e2227..012cf156f628749231ecd47f19962ec5a13c8547 100644 --- a/src/hamiltonian/util/randomdir.h +++ b/src/hamiltonian/util/randomdir.h @@ -31,25 +31,24 @@ * @tparam Enable A helper to distinguish cardinal types from floating point * types and choose a different implementation. */ -template +template struct Rnddir_Helper { - /** Draw a random directions on an nD unit sphere. + /** Draw a random directions on an n-D unit sphere. * * Formulas are from this * [external PDF](https://sites.math.washington.edu/~morrow/335_12/sphericalCoords.pdf) with some hints. * @note For the used rng, we assume it behaves like the RNGCache ... * @param rn where to get the random numbers from. - * @return a random direction uniformly distributed on the unit sphere. + * @param retval On input: storage to use for the new random vector. On return a random direction uniformly distributed on the unit sphere. */ - static auto rnddir(RNG& rn) -> std::array + static void set_to_rnddir(RNG& rn, T& retval) noexcept { - std::array retval; - retval[0] = 1;//beginning of recursion + //beginning of recursion already defined outside. if(SymD > 1) //poor mans constexpr if { //set up that auxiliary data that might be erased later - valuetype angles[SymD-1]; + typename T::value_type angles[SymD-1]; for(int i = 0; i < SymD - 1; ++i) angles[i] = 2.0*rn.real(); //recursion for the polar(?) angles @@ -63,7 +62,6 @@ struct Rnddir_Helper retval[SymD - 1] = retval[SymD - 2] * std::cos(M_PI*angles[SymD - 2]); retval[SymD - 2] = retval[SymD - 2] * std::sin(M_PI*angles[SymD - 2]); } - return retval; } }; @@ -76,20 +74,18 @@ struct Rnddir_Helper */ template struct Rnddir_Helper::value>::type> // enable only for arithmetic types +typename std::enable_if::value>::type> // enable only for arithmetic types { /** Draw a random 1D direction. * * So effectively fot all arithmetic types, this is either +1 or -1. * @note For the used rng, we assume it behaves like the RNGCache ... * @param rn where to get the random numbers from. - * @return Either +1 or -1, chosen randomly. + * @param t where to write the bit. On return: Either +1 or -1, chosen randomly. */ - static auto rnddir(RND& rn) -> std::array + static void set_to_rnddir(RND& rn, T& t) noexcept { - std::array retval = {T(1)}; - if (rn.real() < 0.5) retval[0] = -T(1); - return retval; + if (rn.real() < 0.5) t[0] = -t[0]; } }; @@ -99,11 +95,13 @@ typename std::enable_if::value>::type> // enable only for * * @see Rnddir_Helper * @param rn the RNG to draw the required random numbers from. + * @param t * @return the random direction in the respective space. */ -template -auto rnddir(RNG& rn) -> std::array +template +void set_to_rnddir(RNG&& rn, T& t) noexcept { - return Rnddir_Helper::rnddir(rn); + t[0] = typename T::value_type{1}; + Rnddir_Helper::value>::set_to_rnddir(std::forward(rn), t); } #endif diff --git a/src/libmarqov/core.h b/src/libmarqov/core.h index d33791b2b5d21ba3124994f4c00e015316afb6c2..1bef375a8b4c95e5735a56a365697a6a52f7fadb 100644 --- a/src/libmarqov/core.h +++ b/src/libmarqov/core.h @@ -869,7 +869,8 @@ class Core : public RefType constexpr int SymD = std::tuple_size::value; for (decltype(this->grid.size()) i = 0; i < this->grid.size(); ++i) { - statespace[i] = rnddir, typename StateVector::value_type, SymD>(rngcache); + set_to_rnddir(rngcache, statespace[i]); +// statespace[i] = rnddir, typename StateVector::value_type, SymD>(rngcache); } } diff --git a/src/libmarqov/metropolishelpers.h b/src/libmarqov/metropolishelpers.h index 0695a8ab955726308cb3d5ba39b086145321be08..f71226a6f973f137b20d0d55f3a040dfb576fc8c 100644 --- a/src/libmarqov/metropolishelpers.h +++ b/src/libmarqov/metropolishelpers.h @@ -285,7 +285,7 @@ namespace MARQOV * @param idx the site to be considered * @param hasterms encodes whether the lattice provides the termselector function * - * @return integer sequence (0,1,2,...) with lengt corresponding to number of on-site terms + * @return integer sequence (0,1,2,...) with length corresponding to number of on-site terms */ template std::vector get_terms(const Grid& grid, const Hamiltonian& ham, int idx, std::false_type hasterms) @@ -402,7 +402,7 @@ namespace MARQOV // gather neighbours const auto nbrs = getnbrs(grid, a, rsite); - StateVector neighbourhood = {0}; + StateVector neighbourhood{0}; // sum over neighbours for (std::size_t i=0; i +#include +#include +#include + +#include + +auto dot(const Eigen::Matrix& a, const Eigen::Matrix& b) +{ + return a.dot(b); +} + +//inject eigen and tuple_size compatibility into the standard namespace. +namespace std { + template < + typename _Scalar, + int _Rows, int _Cols, + int _Options, + int _MaxRows, int _MaxCols + > class tuple_size< + Eigen::Matrix< + _Scalar, + _Rows, _Cols, + _Options, + _MaxRows, _MaxCols + > + > { + public: + static_assert( + _Rows != Eigen::Dynamic and _Cols != Eigen::Dynamic, + "tuple_size is only supported for fixed size matrices" + ); + static const size_t value = _Rows * _Cols; + }; +}; + +//include the MARQOV library +#include "libmarqov/libmarqov.h" + +//include the RegularLattice +#include "lattice/regular_hypercubic.h" + +//include some predefined observables, e.g. the magnetization and the energy +#include "hamiltonian/util/observables.h" + + + +// Define interaction term for the Heisenberg model +class MyHeisenberg_interaction +{ + public: + MyHeisenberg_interaction(const double J) : J(J) {} + Eigen::Vector3d get (const Eigen::Vector3d& phi) {return phi;}; + const double J; +}; + +class MyEigenHeisenberg +{ + public: + + // The spin dimension of the Heisenberg or O(3) model + constexpr static int SymD = 3; + + // Define the state vector that this model will use. + // For the Heisenberg model this is a three-dimensional unit vector + typedef Eigen::Vector3d StateVector; + + // Parameters + double J; // The coupling constant + const std::string name; // every Hamiltonian MUST have a name, this is required for the HDF5 output + + // Hamiltonian terms + // here this is only the canonical O(3) interaction, defined above + std::array interactions = {new MyHeisenberg_interaction(J)}; + + // Constructor + MyEigenHeisenberg(double J) : J(J), name("MyEigenHeisenberg"), obs_e(*this){} + ~MyEigenHeisenberg() {delete interactions[0];} + + // Observables + Magnetization obs_m; + Energy obs_e; + decltype(std::make_tuple(obs_m, obs_e)) observables = {std::make_tuple(obs_m, obs_e)}; +}; + + + +using namespace std; +using namespace MARQOV; + +int main() +{ + // Initialize the lattice + int L = 8; + int dim = 2; + RegularHypercubic mylatt(L, dim); + + // Set Monte Carlo parameters using MARQOV::Config + MARQOV::Config mp("./"); // output path + mp.setnmetro(5); // number of Metropolis sweeps per EMCS + mp.setncluster(10); // number of Wolff updates per EMCS + mp.setwarmupsteps(500); // number of EMCS for warmup + mp.setgameloopsteps(3000); // number of EMCS for production + + // Set the Hamiltonian parameters, J, and the inverse temperature beta + double beta = 0.66; + double J = 1; + auto hp = make_tuple(beta, J); + + // Let's create some parameters tuples, a temperature scan for the scheduler to work on + + // Prepare argument list, contains + // 1) reference to the lattice + // 2) Monte Carlo parameter object + // 3) Hamiltonian parameters packed as tuple + auto args = make_tuple(std::ref(mylatt), mp, hp); + + // Eexecute the Core routines of MARQOV. + auto mysim = makeCore(args); + mysim.init(); + mysim.wrmploop(); + mysim.gameloop(); +}