Application Center - Maplesoft

# Groundstate energy of Heisenberg model

You can switch back to the summary page by clicking here.

sample.mws

WorkSheet For Calculating Groundstate Energy of Heisenberg Model

Marvin Weinstein, Stanford Linear Accelerator Center (SLAC)

My goal is to evaluate the CORE approximation formula for the groundstate energy density of a 1+1-dimensional Heisenberg antiferromagnet, based upon a single-state truncation procedure out to range-6 and see how well the basic cluster expansion converges and then, how it can be accelerated using type-II Pade approximants. I will also be evaluating how well I can get a few low lying states in a large set of block from keeping a limited number of states in sub-blocks.

> restart;with(linalg);with(LinearAlgebra);infolevel[LinearAlgebra] := 0:

```Warning, the protected names norm and trace have been redefined and unprotected
```

```Warning, the assigned name GramSchmidt now has a global binding
```

Before beginning the calculations I will create a few routines for taking tensor products of vectors and matrices. This will allow me to build big lattice Hamiltonians in a straightforward manner.

First a general, but somewhat slower technique for taking the tensor product of two arbitrary matrices or vectors. This generic algorithm will be useful however at later stages of the computation.

NB. One could try filling matrices using lists, however this is very slow.

Presumably this has to do with the way Maple 6 exports internal data structures to NAG.

> `&T` := proc(A::Matrix,B::Matrix)
local i,j,nra,nrb,nca,ncb;
nra := RowDimension(A); nrb := RowDimension(B);
nca := ColumnDimension(A); ncb := ColumnDimension(B);
Matrix(1..nra*nrb,1..nca*ncb,proc(i,j) A[iquo(i-1,nrb)+1,iquo(j-1,ncb)+1]*
B[i-iquo(i-1,nrb)*nrb,j-iquo(j-1,ncb)*ncb] end,datatype=float,shape=hermitian);
end:

Next a routine which sorts the eigenvalues returned by Eigenvalues() and gives the smallest one.

Note that eigenvalues are returned in a form which min can't operate on so it is necessary to convert them to a list, make them real and then take the minimum.

> mineig := proc(X) min(op(map(Re,convert(X,list)))) end:

Here is a routine which cleans up printing of things with small numbers resulting from the

fact that Digits = 9 by default

> setnull := proc(x,acc)
if abs(x) < acc then 0
else x
end if
end:

Next, I define the basic building blocks of the one-site problem. The generic Hamiltonian will be

I begin by defining a few basic matrices.

> Sx := 0.5*Matrix(1..2,1..2,[[0.0,1.0],[1.0,0.0]],shape=hermitian);

> Sy := 0.5*Matrix(1..2,1..2,[[0.0,-1.0*I],[1.0*I,0.0]],shape=hermitian);

> Sz := 0.5*Matrix(1..2,1..2,[[1.0,0.0],[0.0,-1.0]],shape=hermitian);

> sxx := Sx &T Sx;

> syy := Matrix(map(Re,Sy &T Sy),shape=hermitian);

> szz := Matrix(map(Re,Sz &T Sz),shape=hermitian);

> umat := IdentityMatrix;

At this point I will define the Hamiltonian for a two-site lattice and find its eigenstates.

> H2 := Matrix(4,4,shape=hermitian);

To speed things up and use as little memory as possible, in order to be able to do a 12-site lattice by brute force, we use the inplace option and define matrices to

be shape=hermitian

> energies2:=Eigenvalues(H2);

> eps2conn := map(Re,mineig(energies2));

Now, to create the Hamiltonian on a big sublattice it is important to observe that it is a sum of tensor products of H2 by unit matrices to the left and right. This special form is particularly quick to construct if one uses special routines which don't bother writing zeros into the Matrix.

The following two routines do this quick tensor product. MakeMatrixUH constructs the operator which is the tensor product of p-unit matrices (2x2) appearing to the left of H2.

Here is where specifying datatype=float and shape=hermitian really speeds things up and saves a lot of memory.

> MakeMatrixUH := proc(p::integer,H::Matrix)
local rH,rU,i,Temp,rTemp;
rU := 2^p;
rH := RowDimension(H);
rTemp := rU*rH;
Temp := Matrix(rTemp,rTemp,datatype=float,shape=hermitian);
for i from 1 to rU do
Temp[(i-1)*rH+1..i*rH,(i-1)*rH+1..i*rH]:=H[1..rH,1..rH]:
od:
RETURN(Temp):
end:

The routine MakeMatrixHU constructs the tensor product of H2 with p-unit matrices appearing to its right. Together MakeMatrixUH and MakeMatrixHU are all we need to construct an arbitrary multi-site Hamiltonian.

> MakeMatrixHU := proc(H::Matrix,p::integer)
local i1,j1,rH,rU,rTemp,Temp;
rH := RowDimension(H);
rU := 2^p;
rTemp := rU*rH;
Temp := Matrix(rTemp,rTemp,datatype=float,shape=hermitian);
for i1 from 1 to rTemp by rU do
for j1 from 1 to rTemp by rU do
Temp[i1..i1+rU-1,j1..j1+rU-1] := H2[iquo(i1-1,rU)+1,iquo(j1-1,rU)+1]*IdentityMatrix(rU);
od;
od;
RETURN(Temp);
end:

Now let us do a three site problem. Note that I call Add explicitly with the inplace option so as to save a lot of memory and time.

> map(setnull,Eigenvalues(H3),1e-10);

> H3;

Even though for the CORE formula we will only need the groundstate energy for 2, 4, 6, 8, 10 and 12-site blocks this exercise is useful to show that things work. We know the answer should be one spin-3/2 multiplet of energy .5, one spin-1/2 multiplet of energy 0 and one spin-1/2 multiplet of energy -1. As we see, the Maple program gave the correct answer. Next let us do the four-site problem and construct along the way the range-2 connected contirbution to the total energy density and the range-2 approximation to the total energy.

> energies4 := map(setnull,Eigenvalues(H4),1e-10);

> en4 := mineig(energies4);

> eps4conn := en4 - 2*eps2conn;

> er1 := eps2conn/2; er2 := (eps2conn+eps4conn)/2;

Next let us construct the 6-site Hamiltonian and do the same computation.

> energies6 := map(setnull,Eigenvalues(H6),1e-10);

> en6 := mineig(energies6);

> eps6conn := en6 - 2*eps4conn - 3*eps2conn;

> er3 := (eps2conn+eps4conn+eps6conn)/2;

> energies8 := map(setnull,map(Re,Eigenvalues(H8)),1e-9);

> en8 := mineig(energies8);

> eps8conn := en8 - 2*eps6conn - 3*eps4conn -4*eps2conn;

> er4 := (eps2conn+eps4conn+eps6conn+eps8conn)/2;

> exact := -.4431472;

> save eps2conn,eps4conn,eps6conn,eps8conn,Sx,Sy,Sz,H2,H4,H5,H6,H8,H10,er1,er2,er3,er4,exact,"heisvals.m";

> save `&T`,MakeMatrixHU,MakeMatrixUH,setnull,mineig,umat,"Mtools.m";

> energies10 :=map(setnull,map(Re,Eigenvalues(H10)),1e-10);

> en10 := mineig(energies10);

> eps10conn := en10-2*eps8conn-3*eps6conn-4*eps4conn-5*eps2conn;

> er5 := (eps2conn+eps4conn+eps6conn+eps8conn+eps10conn)/2;

> exact;

> save eps2conn,eps4conn,eps6conn,eps8conn,eps10conn,Sx,Sy,Sz,H2,H4,H5,H6,H8,H10,er1,er2,er3,er4,er5,exact,"heisvals.m";

The series for the groundstate energy density is quite convergent. As a function of the

range is behaves like ... , where n stands for the range of the

computation. Cleary we are interested in the behavior of this series for large n or

. The best way to take this limit is to fit the computed values for n=1,2,3,...

to a ratio of two polynomials; i.e., for example

etc. One then solves for the coefficients of the rational polynomial which satisfy these equations and extracts A_0 . The program pade2 takes a list of energies and two integers m and n which

specify the maximumum powers of 1/n appearing in the numerator and denominator

respectively.

local A,B,i,r,j,oneoverA,oneoverB,tt,Aops,Bops,eqs,ans;
Aops := a0,seq(cat('a',j),j=2..m+1);
Bops:=seq(cat('b',j),j=2..n+1);
A := [Aops]; B:= [1,Bops];
oneoverA := [1,seq(1/tt^j,j=2..m+1)];
oneoverB := [1,seq(1/tt^j,j=2..n+1)];
ans := solve({eqs},{Aops,Bops});
subs(ans,a0);
end:

>

> epsenv := [er1,er2,er3,er4,er5];

> exact;

Before going any further let's save the functions and intermediate calculations since now we are going to make really big matrices and could run out of memory and crash.

To get the range-6 calculation we need to do the 12-site problem. First we construct the 11-site Hamiltonian and then use it to construct the 12-site one.

The big work is now done, so save the intermediate results. Then, let's get the

eigenvalues of the 12-site Hamiltonian.

> save eps2conn,eps4conn,eps6conn,eps8conn,Sx,Sy,Sz,H2,H4,H5,H6,H8,H10,H11,H12,er1,er2,er3,er4,exact,"heisvals.m";

> e12 := Eigenvalues(H12);

> en12 := mineig(e12);

At this point all that remains is to compute the interesting Pade approximants and

compare the result to the exact answer.

> eps12conn := en12 - 2*eps10conn-3*eps8conn-4*eps6conn-5*eps4conn-6*eps2conn;

> eps10conn;

> er6 := (eps2conn+eps4conn+eps6conn+eps8conn+eps10conn+eps12conn)/2;

> er1;er2;er3;er4;er5;er6;

> exact;

> save eps2conn,eps4conn,eps6conn,eps8conn,eps10conn,eps12conn,Sx,Sy,Sz,H2,H4,H5,H6,H8,H10,H11,H12,er1,er2,er3,er4,er5,er6,exact,"heisvals.m";

> epsenv;

> epsenv := [op(epsenv),er6];