diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2021-11-09 17:59:53 +0100 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2021-11-09 17:59:53 +0100 |
commit | 58e62d50243bc45dbd4ae466642e801a1b46dd69 (patch) | |
tree | 609e60195eff184f43b86d927b2fb6191016837e /collectStokesData.hpp | |
parent | bcc5c327016a3d188a9bdcb0069cc63b0bb8163f (diff) | |
download | code-58e62d50243bc45dbd4ae466642e801a1b46dd69.tar.gz code-58e62d50243bc45dbd4ae466642e801a1b46dd69.tar.bz2 code-58e62d50243bc45dbd4ae466642e801a1b46dd69.zip |
Progress towards a systematized approach.
Diffstat (limited to 'collectStokesData.hpp')
-rw-r--r-- | collectStokesData.hpp | 52 |
1 files changed, 52 insertions, 0 deletions
diff --git a/collectStokesData.hpp b/collectStokesData.hpp new file mode 100644 index 0000000..50aad24 --- /dev/null +++ b/collectStokesData.hpp @@ -0,0 +1,52 @@ +#include "stokes.hpp" +#include <iostream> + +using Complex = std::complex<Real>; + +template<int ...ps, class Generator, typename... T> +void collectStokesData(unsigned N, Generator& r, double ε, T... μs) { + pSpinModel<Real, ps...> ReM(N, r, μs...); + + std::normal_distribution<Real> Red(0, 1); + + Vector<Real> xMin = randomMinimum(ReM, Red, r, ε); + + Real Hx; + Vector<Real> dHx; + Matrix<Real> ddHx; + std::tie(Hx, dHx, ddHx, std::ignore) = ReM.hamGradHess(xMin); + Eigen::EigenSolver<Matrix<Real>> eigenS(ddHx - xMin.dot(dHx) * Matrix<Real>::Identity(xMin.size(), xMin.size()) / Real(N)); + + complex_normal_distribution<Real> d(0, 1, 0); + + pSpinModel M = ReM.template cast<Complex>();; + Vector<Complex> zMin = xMin.cast<Complex>(); + + Vector<Complex> zSaddle = zMin; + + while ((zSaddle - zMin).norm() < 1e-3 * N || abs(imag(M.getHamiltonian(zSaddle))) < 1e-10) { + Vector<Complex> z0 = normalize(zSaddle + 0.5 * randomVector<Complex>(N, d, r)); + zSaddle = findSaddle(M, z0, ε); + } + + Complex H1 = M.getHamiltonian(zMin); + Complex H2 = M.getHamiltonian(zSaddle); + + Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real()); + + M *= exp(Complex(0, φ)); + + Cord c(M, zMin, zSaddle, 3); + c.relaxNewton(10, 1, 1e4); + + ((std::cout << ps), ...) << std::endl;; + ((std::cout << μs), ...) << std::endl;; + + std::apply([](const Tensor<Real, ps>&... Js) -> void { + ((std::cout << Js << std::endl), ...); + } , ReM.Js); + + std::cout << xMin.transpose() << std::endl; + std::cout << Hx << std::endl; + std::cout << eigenS.eigenvalues().real().transpose() << std::endl; +} |