1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
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;
}
|
