examples/ShellMatrix_petsc_Hamiltonian

examples/ShellMatrix_petsc_Hamiltonian/src/main.cpp

  1#include <iostream>
  2#include <mpi.h>
  3
  4/* slepceps includes eigen value routines */
  5#include <slepcmfn.h> 
  6#include <slepceps.h> 
  7
  8/* Maximal system size */
  9#define MaxSites 128
 10
 11/* Hilbert space dimension */
 12#define Q 2
 13
 14/* ScalarType */
 15#define ScalarType double
 16
 17#include "DanceQ.h"
 18using namespace danceq;
 19
 20int main(int argc, char *argv[]) {
 21
 22
 23
 24	PetscInitialize(&argc, &argv, NULL, NULL);
 25
 26	/* System size */
 27	PetscInt L = 14;
 28
 29	/* Particle number */
 30	PetscInt n = 1;
 31
 32	/* Anisotropy */
 33	PetscReal Jz = 0.5;
 34
 35
 36	PetscOptionsGetInt(NULL, NULL, "-L", &L, NULL);
 37	PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
 38	PetscOptionsGetReal(NULL, NULL, "-Jz", &Jz, NULL);
 39
 40	/* Number of MPI ranks */
 41	PetscMPIInt world_size;
 42
 43	/* Rank number */
 44	PetscMPIInt myrank;
 45
 46	PetscCall(SlepcInitialize(&argc, &argv, (char *)0, (char *)0));
 47	PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &world_size));
 48	PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &myrank));
 49
 50	if(myrank == 0){
 51		std::cout << "[myrank: " << myrank << "] - [main] - world_size = " << world_size << std::endl;
 52		std::cout << "[myrank: " << myrank << "] - [main] - L = " << L << std::endl;
 53		std::cout << "[myrank: " << myrank << "] - [main] - n = " << n << std::endl;
 54		std::cout << "[myrank: " << myrank << "] - [main] - Jz = " << Jz << std::endl;
 55	}
 56
 57
 58	/* Operator class */
 59	Hamiltonian_U1 H(L,n); 
 60
 61	/* Simply XXZ-Heisenberg model with periodic boundary conditions */
 62	for(uint64_t i = 0; i < L; i++){
 63		H.add_operator(.5, {i,(i+1)%L}, {"S+","S-"});
 64		H.add_operator(.5, {i,(i+1)%L}, {"S-","S+"});
 65		H.add_operator(Jz, {i,(i+1)%L}, {"Sz","Sz"});
 66	}
 67
 68	uint64_t dim = H.get_dim();
 69	if(myrank == 0){
 70		std::cout << "[myrank: " << myrank << "] - [main] - dim = " << dim << std::endl;
 71		H.info();
 72	}
 73
 74	/* Creating the matrix shell for matrix-free for multiplication */
 75	auto H_matrix = H.create_PetscShellMatrix(true);
 76
 77
 78	/* Setting up the eigenvalue problem */
 79	EPS eps;
 80	PetscCall(EPSCreate(PETSC_COMM_WORLD, &eps));
 81	PetscCall(EPSSetOperators(eps, H_matrix, NULL));
 82	PetscCall(EPSSetProblemType(eps, EPS_HEP));
 83	PetscCall(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL));
 84	PetscCall(EPSSetDimensions(eps, 1, PETSC_DEFAULT, PETSC_DEFAULT));
 85
 86	/* Calculate ground state energy */
 87	PetscCall(EPSSolve(eps));
 88	PetscInt nconv;
 89	PetscCall(EPSGetConverged(eps, &nconv));
 90	PetscScalar E0, E0_imag;
 91	PetscCall(EPSGetEigenvalue(eps,0, &E0, &E0_imag));
 92
 93	if(myrank == 0){
 94		std::cout << "[myrank: " << myrank << "] - [main] - ground state energy = " << E0 << std::endl;
 95	}
 96
 97
 98	PetscCall(MatDestroy(&H_matrix));
 99
100	PetscCall(PetscFinalize());
101
102    return 0;
103}