1 files changed, 11 insertions, 6 deletions
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();