examples/SparseMatrix_petsc_Hamiltonian_imaginary_time_evolution

examples/SparseMatrix_petsc_Hamiltonian_imaginary_time_evolution/src/main.cpp

  1#include <iostream>
  2#include <mpi.h>
  3/* slepcmfn includes matrix exponential */
  4#include <slepc.h> 
  5#include <petsc.h>
  6
  7
  8/* Maximal system size */
  9#define MaxSites 64
 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	/* Initializes Petsc and MPI */
 23	PetscInitialize(&argc, &argv, NULL, NULL);
 24
 25	/* System size */
 26	PetscInt L = 14;
 27
 28	/* Number of steps */
 29	PetscInt steps = 10;
 30
 31	/* Size of step */
 32	PetscReal delta_beta = 0.5;
 33
 34	/* Field strength */
 35	PetscReal h = 0.5;
 36
 37
 38	/* Input using Petsc */
 39	PetscOptionsGetInt(NULL, NULL, "-L", &L, NULL);
 40	PetscOptionsGetInt(NULL, NULL, "-steps", &steps, NULL);
 41	PetscOptionsGetReal(NULL, NULL, "-delta_beta", &delta_beta, NULL);
 42	PetscOptionsGetReal(NULL, NULL, "-h", &h, NULL);
 43
 44
 45	/* Number of MPI ranks */
 46	PetscMPIInt world_size;
 47
 48	/* Rank number */
 49	PetscMPIInt myrank;
 50
 51	PetscCall(SlepcInitialize(&argc, &argv, (char *)0, (char *)0));
 52	PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &world_size));
 53	PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &myrank));
 54
 55	if(myrank == 0){
 56		std::cout << "[myrank: " << myrank << "] - [main] - world_size = " << world_size << std::endl;
 57		std::cout << "[myrank: " << myrank << "] - [main] - delta_beta = " << delta_beta << std::endl;
 58		std::cout << "[myrank: " << myrank << "] - [main] - steps = " << steps << std::endl;
 59		std::cout << "[myrank: " << myrank << "] - [main] - h = " << h << std::endl;
 60	}
 61
 62
 63	/* Hamiltonian */
 64	Hamiltonian_NoU1 H(L); 
 65
 66	/* Defining the transverse field Ising model with J=1 and h */
 67	double J = 1.;
 68	for(uint64_t i = 0; i < L; i++){
 69		H.add_operator(J, {i,(i+1)%L}, {"Sz","Sz"});
 70		H.add_operator(h, {i}, {"Sx"});
 71	}
 72
 73	uint64_t dim = H.get_dim();
 74	if(myrank == 0){
 75		std::cout << "[myrank: " << myrank << "] - [main] - dim = " << dim << std::endl;
 76		H.info();
 77	}
 78
 79	/* Generating the matrices */
 80	auto H_matrix = H.create_PetscSparseMatrix();
 81	if(dim < 10){
 82		PetscCall(MatView(H_matrix, PETSC_VIEWER_STDOUT_WORLD));
 83	}
 84
 85
 86	/* Generating a random initial state */
 87	Vector MyVector(H_matrix);
 88	MyVector.make_random();
 89	Vec Psi = MyVector.create_PetscVec();
 90
 91	/* Second Vec for computation */
 92	Vec PsiPrime;
 93	PetscCall(VecDuplicate(Psi, &PsiPrime));
 94
 95	/* Third Vec for the ground state */
 96	Vec Psi0;
 97	PetscCall(VecDuplicate(Psi, &Psi0));
 98
 99
100
101	/* Computation of the ground state */
102	EPS eps;
103	PetscCall(EPSCreate(PETSC_COMM_WORLD, &eps));
104	PetscCall(EPSSetOperators(eps, H_matrix, NULL));
105	PetscCall(EPSSetProblemType(eps, EPS_HEP));
106	PetscCall(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL));
107	PetscCall(EPSSetDimensions(eps, 1, PETSC_DEFAULT, PETSC_DEFAULT));
108
109	/* Retrieving ground state energy */
110	PetscCall(EPSSolve(eps));
111	PetscInt nconv;
112	PetscCall(EPSGetConverged(eps, &nconv));
113	if(nconv == 0){
114		if(myrank == 0){
115			std::cout << "[ERROR] - [myrank: " << myrank << "] - [main] - ground state not converged" << std::endl;
116		}
117		return 1;
118	}
119	PetscScalar E0;
120	PetscCall(EPSGetEigenpair(eps,0, &E0, NULL, Psi0, NULL));
121
122	if(myrank == 0){
123		std::cout << "[myrank: " << myrank << "] - [main] - ground state energy = " << E0 << std::endl;
124	}
125
126
127
128	/* Shifting the spectrum of H such that all eigenvalues are semi positive */
129	H.add_diagonal(-E0);
130
131	/* Extracting the new Hamiltonian */
132	H_matrix = H.create_PetscSparseMatrix();
133
134
135	/* Setting up the matrix exponential with delta_beta */
136	MFN mfn;
137	FN f;
138	PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
139	PetscCall(MFNSetOperator(mfn,H_matrix));
140	PetscCall(MFNGetFN(mfn,&f));
141	PetscCall(FNSetType(f,FNEXP));
142	PetscCall(FNSetScale(f,-delta_beta,1.0)); /* Setting the delta_beta steps */
143	PetscCall(MFNSetTolerances(mfn,1e-07,PETSC_DEFAULT));
144	PetscCall(MFNSetFromOptions(mfn));
145
146
147	/* Loop of steps */
148	double beta = 0.;
149	PetscScalar energy,overlap;
150	PetscReal norm;
151	for(uint64_t i = 0; i < steps; i++){
152
153		beta += delta_beta;
154
155		/* Time evolution PsiPrime = e^{-delta_beta (H-E0)} Psi */
156		PetscCall(MFNSolve(mfn,Psi,PsiPrime));
157
158		/* Setting Psi = PsiPrime */
159		PetscCall(VecCopy(PsiPrime, Psi)); 
160
161		/* Calculating the norm */
162		PetscCall(VecNorm(Psi,NORM_2,&norm));
163
164		/* Rescale Psi */
165		PetscCall(VecScale(Psi, 1./norm));
166
167		/* Calculating the overlap with the ground state */
168		PetscCall(VecDot(Psi, Psi0, &overlap));
169
170		/* Perform matrix-vector multiplication: PsiPrime = (H-E0) * Psi */
171		PetscCall(MatMult(H_matrix, Psi, PsiPrime));
172
173		/* Calculating the energy */
174		PetscCall(VecDot(Psi, PsiPrime, &energy));
175
176		/* Output */
177		if(myrank == 0){
178			std::cout << "[myrank: " << myrank << "] - [main] - beta = " << beta << ": |<Psi|Psi0>| = " << std::abs(overlap) << " - <Psi|H-E0|Psi> = " << energy << std::endl;
179		}
180
181	}
182
183	/* Deleting Petsc objects */
184	PetscCall(MatDestroy(&H_matrix));
185	PetscCall(VecDestroy(&Psi));
186	PetscCall(VecDestroy(&Psi0));
187	PetscCall(VecDestroy(&PsiPrime));
188
189	/* Cleaning up */
190	PetscCall(PetscFinalize());
191
192    return 0;
193}