diff options
| -rw-r--r-- | langevin.cpp | 16 |
1 files changed, 9 insertions, 7 deletions
diff --git a/langevin.cpp b/langevin.cpp index d704266..119a64b 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -13,13 +13,13 @@ #include "complex_normal.hpp" #include "tensor.hpp" -#define P_SPIN_P 3 -const unsigned p = P_SPIN_P; +#define PSPIN_P 3 +const unsigned p = PSPIN_P; // polynomial degree of Hamiltonian using Scalar = std::complex<double>; using Vector = Eigen::VectorXcd; using Matrix = Eigen::MatrixXcd; -using Tensor = Eigen::Tensor<Scalar, p>; +using Tensor = Eigen::Tensor<Scalar, PSPIN_P>; using Rng = randutils::random_generator<pcg32>; @@ -40,9 +40,11 @@ std::tuple<Scalar, Vector, Matrix> hamGradHess(const Tensor& J, const Vector& z) Matrix Jz = contractDown(J, z); Vector Jzz = Jz * z; - Matrix hessian = ((p - 1) * (double)p / factorial(p)) * Jz; - Vector gradient = ((double)p / factorial(p)) * Jzz; - Scalar hamiltonian = (1.0 / factorial(p)) * Jzz.dot(z); + double f = factorial(p); + + Matrix hessian = ((p - 1) * p / f) * Jz; + Vector gradient = (p / f) * Jzz; + Scalar hamiltonian = (1 / f) * Jzz.dot(z); return {hamiltonian, gradient, hessian}; } @@ -142,7 +144,7 @@ int main(int argc, char* argv[]) { Rng r; complex_normal_distribution<> d(0, σ, κ); - Tensor J = generateCouplings<Scalar, p>(N, d, r.engine()); + Tensor J = generateCouplings<Scalar, PSPIN_P>(N, d, r.engine()); Vector z = initializeVector(N, r); std::function<bool(double, unsigned)> findSaddle = [δ](double W, unsigned) { |