diff options
| -rw-r--r-- | collectStokesData.hpp | 52 | ||||
| -rw-r--r-- | dynamics.hpp | 4 | ||||
| -rw-r--r-- | langevin.cpp | 17 |
3 files changed, 65 insertions, 8 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; +} diff --git a/dynamics.hpp b/dynamics.hpp index 36b2dfa..4f33c1a 100644 --- a/dynamics.hpp +++ b/dynamics.hpp @@ -9,7 +9,7 @@ template <class Real, int ...ps> Vector<Real> findMinimum(const pSpinModel<Real, ps...>& M, const Vector<Real>& z0, Real ε) { Vector<Real> z = z0; - Real λ = 10; + Real λ = 100; Real H; Vector<Real> dH; @@ -91,7 +91,7 @@ Vector<Real> randomMinimum(const pSpinModel<Real, ps...>& M, Distribution d, Gen bool foundSaddle = false; while (!foundSaddle) { - Vector<Real> z0 = normalize(randomVector<Real>(M.dimension(), d, r.engine())); + Vector<Real> z0 = normalize(randomVector<Real>(M.dimension(), d, r)); try { zSaddle = findMinimum(M, z0, ε); diff --git a/langevin.cpp b/langevin.cpp index 0e5a616..ea28d65 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -11,6 +11,8 @@ #include "dynamics.hpp" #include "stokes.hpp" +#include "collectStokesData.hpp" + #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" #include "unsupported/Eigen/CXX11/src/Tensor/TensorFFT.h" @@ -83,20 +85,23 @@ int main(int argc, char* argv[]) { Rng r; -// pSpinModel<Real, 2, 4> pSpin({J2, J4}); - pSpinModel<Real, 2, 4>pSpin(N, r.engine(), 1, 0.1); + collectStokesData<2, 4>(N, r.engine(), 1e-15, 1.0, 0.01); + + pSpinModel<Real, 2, 4>pSpin(N, r.engine(), 1, 0.01); std::normal_distribution<Real> Red(0, 1); - RealVector zMin = randomMinimum(pSpin, Red, r, ε); + RealVector zMin = randomMinimum(pSpin, Red, r.engine(), ε); Real Hr; RealVector dHr; RealMatrix ddHr; std::tie(Hr, dHr, ddHr, std::ignore) = pSpin.hamGradHess(zMin); - Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size() + zMin * zMin.transpose() / (Real)N); + RealMatrix M1 = ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size(); + RealMatrix M2 = ddHr - zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) / Real(N); + Eigen::EigenSolver<Matrix<Real>> eigenS(M2); for (unsigned i = 0; i < N; i++) { - RealVector zNew = normalize(zMin + 1e-3 * eigenS.eigenvectors().col(i).real()); - std::cout << eigenS.eigenvalues()(i) << " " << 2.0 * (pSpin.getHamiltonian(zNew) - Hr) / 1e-6 << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << " " << eigenS.eigenvectors().col(0).dot(eigenS.eigenvectors().col(i)) << std::endl; + RealVector zNew = normalize(zMin + 1e-5 * eigenS.eigenvectors().col(i).real()); + std::cout << eigenS.eigenvalues()(i) << " " << 2.0 * (pSpin.getHamiltonian(zNew) - Hr) / 1e-10 << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << " " << eigenS.eigenvectors().col(0).dot(eigenS.eigenvectors().col(i)) << std::endl; } std::cout << std::endl; getchar(); |