examples/ShellMatrix

examples/ShellMatrix/src/main.cpp

  1#include <iostream>
  2#ifdef USE_MPI
  3#include <mpi.h>
  4#endif
  5#ifdef USE_OMP
  6#include <omp.h>
  7#endif
  8
  9/* Including Eigen */
 10#include <Eigen/Dense>
 11#include <Eigen/Eigenvalues>
 12
 13/* Maximal system size */
 14#define MaxSites 128
 15
 16/* Hilbert space dimension */
 17#define Q 2
 18
 19/* ScalarType */
 20#define ScalarType std::complex<double>
 21
 22/* You have to include Eigen before including the DanceQ */
 23#include "DanceQ.h"
 24using namespace danceq;
 25
 26
 27int main(int argc, char *argv[]) {
 28
 29	/* Rank number */
 30	int myrank = 0;
 31#ifdef USE_MPI
 32	/* Initializes MPI */
 33	MPI_Init(&argc, &argv);
 34
 35	/* Number of MPI ranks */
 36	int world_size;
 37	MPI_Comm_size(MPI_COMM_WORLD, &world_size);
 38
 39	/* Rank number */
 40    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
 41#endif
 42
 43	/* Number of particles */
 44	uint64_t n = 1;
 45
 46	/* System size */
 47	uint64_t L = 4;
 48
 49	/* number of states */
 50	uint64_t number_of_states = 4;
 51
 52	/* System size */
 53	double Jz = 4.;
 54
 55	/* Particle number and system size can be edited at run time: mpirun -n NUM_CORES ./main -n 5 -L 10 */
 56	if (argc > 1){
 57		for(uint32_t i = 1; i < argc; ++i){ 
 58			std::string arg_val(argv[i]);
 59			bool input = true;
 60			if(arg_val == "-n" and i < argc - 1){
 61				n = static_cast<uint64_t>(std::stoi(argv[i+1]));
 62				input = false;i++;
 63			}
 64			if(arg_val == "-L" and i < argc - 1){
 65				L = static_cast<uint64_t>(std::stoi(argv[i+1]));
 66				input = false;i++;
 67			}
 68			if(arg_val == "-number_of_states" and i < argc - 1){
 69				number_of_states = static_cast<uint64_t>(std::stoi(argv[i+1]));
 70				input = false;i++;
 71			}
 72			if(arg_val == "-Jz" and i < argc - 1){
 73				Jz = static_cast<double>(std::stof(argv[i+1]));
 74				input = false;i++;
 75			}
 76			if(input and myrank == 0){
 77				std::cout << "[myrank: " << myrank << "] - [main] - Unknown option at input variable " << i << " : " << arg_val << std::endl;
 78			}
 79		}
 80	}
 81
 82
 83	if(myrank == 0){
 84		std::cout << "[main] - L = " << L << std::endl;
 85		std::cout << "[main] - n = " << n << std::endl;
 86		std::cout << "[main] - number_of_states = " << number_of_states << std::endl;
 87		std::cout << "[main] - Jz = " << Jz << std::endl;
 88	}
 89
 90	/* Hamiltonian */
 91	Hamiltonian_U1 H(L,n); 
 92
 93
 94	/* XXZ-Heisenberg model with open boundary conditions and a tilted field */
 95	for(uint64_t i = 0; i < L-1; i++){
 96		H.add_operator(.5, {i,(i+1)%L}, {"S+","S-"});
 97		H.add_operator(.5, {i,(i+1)%L}, {"S-","S+"});
 98		H.add_operator(Jz, {i,(i+1)%L}, {"Sz","Sz"});
 99		H.add_operator(static_cast<double>(i), {i}, {"Sz"});
100	}
101
102	/* Information */
103	H.info();
104
105	/* Observable */
106	Hamiltonian_U1 O(L,n); 
107
108	/* Operator measuring the magnetization of the first site */
109	O.add_operator(1., {0}, {"Sz"});
110
111
112
113	auto O_matrix = O.create_ShellMatrix();
114
115	std::vector<Vector> states;
116	auto data = lanczos(H, number_of_states, &states /* returns eigen states if this is not nullptr */);
117	for(uint64_t s_for = 0UL; s_for < data.size(); s_for++){
118		auto obs = O_matrix.get_expectation_value(states[s_for]);
119		if(myrank == 0){
120			std::cout << "[main] - E_" << s_for << " = " << data[s_for] << ", O_" << s_for << " = " << obs << std::endl;
121		}
122	}
123
124#ifdef USE_MPI
125    MPI_Finalize();
126#endif
127
128
129    return 0;
130}