1#include <iostream>
2#include <complex>
3
4/* Including Eigen */
5#include <Eigen/Dense>
6#include <Eigen/Eigenvalues>
7
8/* Maximal system size */
9#define MaxSites 256
10
11/* Hilbert space dimension */
12#define Q 2
13
14/* ScalarType */
15#define ScalarType std::complex<double>
16
17
18/* You have to include Eigen before including the DanceQ */
19#include "DanceQ.h"
20using namespace danceq;
21
22int main(int argc, char *argv[]) {
23
24 /* Number of particles */
25 uint64_t n = 1;
26
27 /* System size */
28 uint64_t L = 4;
29
30 /* number of subsystems */
31 uint64_t Npart = 3;
32
33 /* Particle number, system size and number of subsystem can be edited at run time: ./main -n 5 -L 10 -Npart 4 */
34 if (argc > 1){
35 for(uint32_t i = 1; i < argc; ++i){
36 std::string arg_val(argv[i]);
37 bool input = true;
38 if(arg_val == "-n" and i < argc - 1){
39 n = static_cast<uint64_t>(std::stoi(argv[i+1]));
40 input = false;i++;
41 }
42 if(arg_val == "-L" and i < argc - 1){
43 L = static_cast<uint64_t>(std::stoi(argv[i+1]));
44 input = false;i++;
45 }
46 if(input){
47 std::cout << "[main] - Unknown option at input variable " << i << " : " << arg_val << std::endl;
48 }
49 }
50 }
51
52 std::cout << "[main] - L = " << L << std::endl;
53 std::cout << "[main] - n = " << n << std::endl;
54
55 /* Operator class */
56 Lindbladian_U1 Lind(L,n);
57
58 uint64_t dim = Lind.get_dim();
59 std::cout << "[main] - dim = " << dim << std::endl;
60
61 /* Simply XXX-Heisenberg for the unitary dynamics */
62 for(uint64_t i = 0; i < L; i++){
63 Lind.add_operator(.5, {i,(i+1)%L}, {"S+","S-"});
64 Lind.add_operator(.5, {i,(i+1)%L}, {"S-","S+"});
65 Lind.add_operator(1., {i,(i+1)%L}, {"Sz","Sz"});
66 Lind.add_operator(std::pow(-1.,i), {i}, {"Sz"});
67 }
68
69 /* Jump operator acting on the first site */
70 Lind.add_jump_operator(1., {0}, {"Sz"});
71
72 /* Printing information of the Hamiltonian */
73 Lind.info();
74
75 std::cout << "[main] - generating the matrix ... " << std::endl;
76 auto Lind_matrix = Lind.create_EigenDenseMatrix();
77
78#ifdef EIGEN_USE_MKL_ALL
79 std::cout << "[main] - using MKL" << std::endl;
80#else
81 std::cout << "[main] - not using MKL" << std::endl;
82#endif
83
84 std::cout << "[main] - calculating the entire spectrum ... " << std::endl;
85
86 /* Non-hermitian complex solver */
87 Eigen::ComplexEigenSolver<decltype(Lind_matrix)> solver;
88
89 /* Computation without eigenvectors */
90 solver.compute(Lind_matrix, false);
91
92 auto eigenvalues = solver.eigenvalues();
93
94 std::cout << "[main] - eigenvalues = " << eigenvalues.transpose() << std::endl;
95 std::cout << "[main] - steady state eigenvalue = " << eigenvalues(0) << std::endl;
96 if(Lind.get_dim() > 1){
97 std::cout << "[main] - decay rate = " << eigenvalues(1) << std::endl;
98 }
99
100 return 0;
101}