Infinite DMRG

This notebook takes you through the writing of an infinite DMRG code for the transverse field Ising Model

Load routines for working with Named Tensors

SetDirectory[NotebookDirectory[]]

NotebookEvaluate[FileNameJoin[{NotebookDirectory[],"MinimalNamedTensors.nb"}]]

Big Picture

This is an iterative algorithm.  We will implement a function which does 1 iteration, and then nest it in some sort of control structure.

We will call the function iterateIDMRG

iterateIDMRG::usage=
"IterateIDMRG runs one step of the Infinite DMRG Algorithm

input: list of rules with at least the following keys: 
\"ZL\" -- a named tensor which gives the action of Z on the rightmost spin of the left environment
\"HL\" -- a named tensor which give H on the left environment 
\"J\", \"h\" -- Hamiltonian Parameters
\"E\" -- estimate of the energy per site 

the following rules are optional:
\"tolerance\" -- truncate singular values below this
\"maxbond\" -- keep at most this size of bond dimension

output: a list containing all of the elements of the input, but updated.  Also
\"L\" -- last left normalized Matrix in MPS
\"singularvalues\" -- list of singular values";

Sample Matrices

To test and debug, we will need some sample input -- these can also be used for the initial iteration

The first environment is just “empty”, consisting of both ZL and HL is just a 1 x1 matrix that is zero -- This has the negative effect of not breaking the up-down symmetry, which means that the  ground state in the ordered phase is a “Schrodinger Cat” -- and requiring a higher bond dimension to describe the state

ZL1={"tensor"->{{0}},"names"->{"t","t'"}};
HL1={"tensor"->{{0}},"names"->{"t","t'"}};

The second environment is that the environment spins are all in the +z direction:  That corresponds to a bond dimension 1 with ZL=1 and HL=0 matrix we can use for both ZL and HL is just a 1x1 matrix that is zero.  This breaks the symmetry, and is  a more efficient starting point

ZL2={"tensor"->{{1}},"names"->{"t","t'"}};
HL2={"tensor"->{{0}},"names"->{"t","t'"}};

The third is that the environment consists of a single spin, which can be up or down -- ZL  = {{1,0},{0,-1}} and HL={{0,h},{h,0}} -- this has the same degeneracy issues as the first

ZL3={"tensor"->pauliZ,"names"->{"t","t'"}};
HL3[h_]={"tensor"->SparseArray[{{0,-h},{-h,0}}],"names"->{"t","t'"}};

Constructing the effective H

Code

Make a function which given ZL, HL, J, and h, constructs the Heff.  After each section, stop and add a Return statement to stop the evaluation, and test if it is working

MakeHeff[ZLin_List,HLin_List,J_,h_]:=
Module[
{
(*Put local variables here *)
ZL,HL,ZR,HR,
Z1,Z2,X1,X2,
IL,I1,I2,IR,
HLxIxIxI,ZxZxIxI,IxZxZxI,IxIxZxZ,IxIxIxHR,
IxXxIxI,IxIxXxI,heff
, len , hefftensor

},
(* For robustness, set the names of the indices in
 ZL and HL to some canonical form, and generate
ZR and HR symmetrically *)
(*Here is how I would do it for the Z's -- The wrapping in SparseArray converts dense into sparse -- which will add efficiencies later*)
ZL={"tensor"->SparseArray[("tensor"/.ZLin)],"names"->{"t","t'"}};
ZR={"tensor"->SparseArray[("tensor"/.ZLin)],"names"->{"s","s'"}};
Return[ZL];  (* As you write code, move this break point, and change what it is returning *)
(* Do the same with HL and HR *)


(*Generate the four matrices Z1, Z2, X1, X2, corresponding to X and Z on the two 
central spins*)
Z1={"tensor"->pauliZ,"names"->{"σ","σ'"}};
Z2={"tensor"->pauliZ,"names"->{"τ","τ'"}};
(* Do the same with X1 and X2 *)


(*Generate the four identity matrices, corresponding to the indices of the 
leftt, spin1, spin2, right -- call them IL, I1, I2, IR 
-- I've started you out by extracting the size of the environment, and doing one of them *)

len=Length["tensor"/.ZLin];
IL=id["t","t'",len];

(* generate the other 3 *)



(*By using Combine and CombineIndices generate the seven tensors, and rewrite them as matrices (ie. combine all the non-primed indices, and all the primed indices)
HL*I*I*I, Z*Z*I*I, I*Z*Z*I,I*I*Z*Z, I*I*I*HR,I*X*I*I,I*I*X*I  *)

HLxIxIxI=Combine[HL,I1,I2,IR];
(* Do the rest *)


(*Extract the tensors from each of these and add them, with appropriate factors of J and h  -- I've done the first couple terms*)
heff=("tensor"/.HLxIxIxI);
heff += ("tensor"/.IxIxIxHR);
heff += -J*("tensor"/.ZxZxIxI);

heff += -J*("tensor"/.ZxZxIxI);
heff += -J*("tensor"/.IxZxZxI);
heff += -J*("tensor"/.IxIxZxZ);
heff += -h*("tensor"/.IxXxIxI);
heff += -h*("tensor"/.IxIxXxI);

(* Turn into a NamedTensor *)
hefftensor={"tensor"->heff,"names"->{"t","t'","σ","σ'","τ","τ'","s","s'"}};
hefftensor=CombineIndices[hefftensor,{"t","σ","τ","s"}];
hefftensor=CombineIndices[hefftensor,{"t'","σ'","τ'","s'"}];

(*Send along the info about the dimensions of the different spaces for when we want to break them apart, and the meaning of the indices.  [If nothing else, handy for debugging -- can remove later if they bother you.] *)
Join[hefftensor,{"subdims"->{len,2,2,len},
"subnames"->{{"t","σ","τ","s"},{"t'","σ'","τ'","s'"}},"h"->h,"J"->J}]
]

Debug/Test

Debug each section of the function before going on -- it is easier to debug as you go -- just keep moving the “Return” Function

Some Space for debugging:

Heff1=MakeHeff[ZL1,HL1,1.2,0.3]

Heff2=MakeHeff[ZL2,HL2,1.2,0.3]

Heff3=MakeHeff[ZL3,HL3[0.3],1.2,0.3]

My results to compare with:

Heff1=MakeHeff[ZL1,HL1,1.2,0.3]

MatrixForm["tensor"/.Heff1]

Eigenvalues[Normal["tensor"/.Heff1]]

Heff2=MakeHeff[ZL2,HL2,1.2,0.3]

MatrixForm["tensor"/.Heff2]

Eigenvalues[Normal["tensor"/.Heff2]]

Depending on the function, you can sometimes make symbolic checks.  This is one of those cases

symbolicHeff2=MakeHeff[ZL2,HL2,J,h]

MatrixForm["tensor"/.symbolicHeff2]

Heff3=MakeHeff[ZL3,HL3[0.3],1.2,0.3]

Sort@Eigenvalues[Normal["tensor"/.Heff3]]

Note the near-degeneracy -- this is a consequence of the Schrodinger-Cat structure -- Causes all sorts of numerical issues -- stick with boundaries 2.

Finding the optimal “Q” matrix

Code

FindQ[Heff_,Energypersiteguess_]:=
Module[
{
(*Local Variables *)
Q,EnergyPerSite,eigsys

},
(* Call eig  -- define Q as the output.  Define Energypersite as half the eigenvalue. *)

Return[Heff];  (* As you write code, move this break point, and change what it is returning *)



(*Return Q  and EnergyPerSite  *)
Join[Q,{"EnergyPerSite"->EnergyPerSite,"subnames"->First[("subnames"/.Heff)],
"subdims"->("subdims"/.Heff),"guess"->Energypersiteguess*2
}]
]

Debug / Test

Debug before moving on. -- Use the trick of adding breakpoints with Return

My results, for comparison:

Q1=FindQ[Heff1,-5]//Chop

Q2=FindQ[Heff2,-5]//Chop

Q3=FindQ[Heff3,-5]//Chop

Get the left tensor and singular values from the Q matrix

I’ll write this one -- since we only have so much class time -- feel free to write your own and compare

decompose[Q_,opts___]:=
Module[{QasTensor,Qindex,len,leftdim,tdim,σdim,
QasMatrix,sv,L},
Qindex=First["names"/.Q];
len=Length["tensor"/.Q];

(* I have hard coded the names of the indices, and the
dimension of the physical Hilbert space --
could more generally pull this from Q  *)

leftdim=Sqrt[len]; (*Dimension of combines "t, σ" index*)
tdim=leftdim/2;
σdim=2;

(* Unpack the indices to make matrix *)
QasMatrix=BreakIndices[Q,Qindex,{"left","right"},{leftdim,leftdim}];

(* Run SVD *)
sv=svd[QasMatrix,{"b1","b2"},opts];
(* Convert Left isometry to original variables *)

L=BreakIndices[sv[[1]],"left",{"t","σ"},{tdim,σdim}];

(* Assemble Output *)
{"L"->L,"SingularValues"->Tr["tensor"/.sv[[2]],List]}]

Verify that you get my results before going on -- it should work, since you are using my routine -- but it is always good to debug

d1=decompose[Q1]

d2=decompose[Q2]

d3=decompose[Q3]

My results

d1=decompose[Q1]

d1b=decompose[Q1,"trunc"->1]

d1c=decompose[Q1,"tolerance"->0.6]

d2=decompose[Q2]

Note how much nicer the entanglement structure is when we have polarized boundaries (Q2 has polarized boundaries, Q1 does not)

d3=decompose[Q3]

Expectation Values

It will be useful to calculate <X> <Z> <ZZ> on the central sites, and the correlation length (in the MPS you would get by infinitely repeating (L1) (L2)) -- you already did these in class

singlesiteexpect[L_,op_,singularvalues_]:=
Module[{(*local variables*)
Lg={"tensor"->("tensor"/.L),"names"->{"t","σ","b1"}},
Lconj={"tensor"->Conjugate[("tensor"/.L)],"names"->{"t","σ'","b2"}},
Lambdasq={"tensor"->SparseArray[{Band[{1,1}]->singularvalues^2}],"names"->{"b1","b2"}},
opg={"tensor"->("tensor"/.op),"names"->{"σ","σ'"}},result},
result=autoContract[Lg,Lconj,opg,Lambdasq];
"tensor"/.result]

Because we are constructing the Schrodinger cat, the first boundary condition gives zero for the expectation value

singlesiteexpect["L"/.d1,{"tensor"->pauliZ},"SingularValues"/.d1]

singlesiteexpect["L"/.d2,{"tensor"->pauliZ},"SingularValues"/.d1]

singlesiteexpect["L"/.d2,{"tensor"->pauliZ},"SingularValues"/.d1]

singlesiteexpect["L"/.d1,{"tensor"->pauliZ},"SingularValues"/.d1]

singlesiteexpect["L"/.d2,{"tensor"->pauliZ},"SingularValues"/.d1]

singlesiteexpect["L"/.d2,{"tensor"->pauliX},"SingularValues"/.d1]

twositeexpect[L_,op1_,op2_,singularvalues_]:=
Module[{(*local variables*)
L1={"tensor"->("tensor"/.L),"names"->{"t","σ1","b1"}},
Lconj1={"tensor"->Conjugate[("tensor"/.L)],"names"->{"t","σ1'","b2"}},
R2={"tensor"->("tensor"/.L),"names"->{"s","σ2","c1"}},
Rconj2={"tensor"->Conjugate[("tensor"/.L)],"names"->{"s","σ2'","c2"}},
Lambda={"tensor"->SparseArray[{Band[{1,1}]->singularvalues}],"names"->{"b1","c1"}},
Lambdaconj={"tensor"->SparseArray[{Band[{1,1}]->singularvalues}],"names"->{"b2","c2"}},
op1g={"tensor"->("tensor"/.op1),"names"->{"σ1","σ1'"}},
op2g={"tensor"->("tensor"/.op2),"names"->{"σ2","σ2'"}}
,result},
result=autoContract[L1,Lconj1,op1g,Lambda,Lambdaconj,R2,Rconj2,op2g];
"tensor"/.result
]

ww=twositeexpect["L"/.d1,{"tensor"->pauliZ},{"tensor"->pauliZ},"SingularValues"/.d1]

ww=twositeexpect["L"/.d1,{"tensor"->pauliZ},{"tensor"->pauliZ},"SingularValues"/.d1]

Construct New HL and New ZL

Code

I’ll do makeZL -- and let you do makeHL.  As before, I will assume I know what indices L will come in with: “t”,”σ”, and “b1”.  If you can’t be guaranteed of their names, then you need a way to pass that information (or be sure there is a canonical order you can use to grab them).  For Hermitian matrices, like HL or ZL, the order of the indices do not matter, so you can freely change the names with  “ZLNew={“tensor”/.ZLOld,”names”->{“newname1”,”newname2”}}”

makeZL[L_]:=
Module[
{
(*local variables*)
Lconj={"tensor"->Conjugate[("tensor"/.L)],"names"->{"t","σ'","b2"}},
Z={"tensor"->pauliZ,"names"->{"σ","σ'"}},ZL
},
ZL=autoContract[L,Lconj,Z];
ZL
]

makeHL[L_,oldHL_,oldZL_,h_,J_,EnergyPerSite_]:=
Module[
{
(*local variables*)

},

]

Debug / Test

Debug before moving on.

(For my debugging I outputted a bunch of the intermediate steps -- you don’t have to)

Boundary Conditions 1

newHL1=makeHL["L"/.d1,(*old HL*) HL1,(* old Z *) ZL1, (*h*) 0.3,(*J*) 1.2,"EnergyPerSite"/.Q1]

newHL1=makeHL["L"/.d1,(*old HL*) HL1,(* old Z *) ZL1, (*h*) 0.3,(*J*) 1.2,"EnergyPerSite"/.Q1]

newZL1=makeZL["L"/.d1]

newZL1=makeZL["L"/.d1]

MatrixForm[(1/2){{1,1},{1,-1}}.("tensor"/.newZL1).{{1,1},{1,-1}} //Chop]

MatrixForm["tensor"/.ZL3]

MatrixForm[(1/2){{1,1},{1,-1}}.("tensor"/.newHL1).{{1,1},{1,-1}} //Chop]

MatrixForm["tensor"/.HL3[0.3]]

Boundary Conditions 2

newHL2=makeHL["L"/.d2,(*old HL*) HL2,(* old Z *) ZL2, (*h*) 0.3,(*J*) 1.2,"EnergyPerSite"/.Q2]

newHL2=makeHL["L"/.d2,(*old HL*) HL2,(* old Z *) ZL2, (*h*) 0.3,(*J*) 1.2,"EnergyPerSite"/.Q2]

newZL2=makeZL["L"/.d2]

newZL2=makeZL["L"/.d2]

"L"/.d2

Boundary Conditions 3

newHL3=makeHL["L"/.d3,(*old HL*) HL3[0.3],(* old Z *) ZL3, (*h*) 0.3,(*J*) 1.2,"EnergyPerSite"/.Q2]

newZL3=makeZL["L"/.d3]

newHL3=makeHL["L"/.d3,(*old HL*) HL3[0.3],(* old Z *) ZL3, (*h*) 0.3,(*J*) 1.2,"EnergyPerSite"/.Q2]

newZL3=makeZL["L"/.d3]

Correlation Length

Given the matrix Q, we can extract a correlation length by imagine making a iMPS which simply repeats Q.  We first find the transfer matrix, find its eigenvalues, and extract the correlation length as the inverse of the Log of the ratio of the two largest eigenvalues. -- I just went ahead and wrote it.  Note, both Q1 and Q2 have a zero correlation length and will presumably cause an error.

corlength[Q_List]:=
Module[{Qmat,Qmatstar,transfer,len=(Length["tensor"/.Q]),dims,origname=First["names"/.Q],eigs},
dims={Sqrt[len/4],4,Sqrt[len/4]};
Qmat=BreakIndices[Q,origname,{"l1","sigs","r1"},dims];
Qmatstar={"tensor"->("tensor"/.Qmat),"names"->{"l2","sigs","r2"}};
transfer=autoContract[Qmat,Qmatstar];
transfer=leftCombineIndices[transfer,{"l1","l2"},"l"];
transfer=rightCombineIndices[transfer,{"r1","r2"},"r"];
eigs=Sort@Eigenvalues["tensor"/.transfer];
1/Log[eigs[[-1]]/eigs[[-2]]]
]

corlength[Q3]

Putting it together

I went ahead and drafted something, but depending on the output of your functions, it may need to be tweaked.  It has a bunch of bells and whistles that you don’t need -- for example it can take an optional “verbose” option where it prints to the screen while it is doing things, and it records the timing of the different pieces.  Useful for profiling -- though as I said before, you need to be realistic about how fast you can make this run.

iterateIDMRG[input_List]:=
Module[
{
(******************** Input ***********************)
ZL="ZL"/.input
(*a named tensor which gives the action of Z on the rightmost spin of the left environment*),
HL="HL"/.input  
(*a named tensor which give H on the left environment *),
J="J"/.input, 
h="h"/.input, 
EnergyPerSite="E" /.input,
tolerance="tolerance"/.input (* truncate singular values below this*),
maxbond="maxbond" /.input (* keep at most this size of bond dimension*),
iteration="iteration"/.input //.("iteration"->0),
shift,
verbose="verbose"/.input //.("verbose"->False),
(******************* Working Variables *************)
Heff,eigres,Q,newEnergyPerSite,decomp,opts,L,SingularValues,
newHL,newZL,message1,message2,message3,time,time2,makeHtime,findQtime,SVDtime,Z,X,ZZ,alltime,time0,HZLtime
},

time0=AbsoluteTime[];
(* Fudge factor for eigensolver -- to make sure that we are below spectrum -- might not need it, but is a good knob to have*)
shift="shift"/.input //."shift"->-Max[J,h]; 
(* Make Heff *)
If[verbose,
message1=PrintTemporary["Constructing H"];
];
time=AbsoluteTime[];
Heff=MakeHeff[ZL,HL,J,h];
time2=AbsoluteTime[];
makeHtime=time2-time;
If[verbose,
NotebookDelete[message1];
message1=PrintTemporary["Constructing H took time: "<>ToString[makeHtime]<>"s"];
];
(* Solve Eigenproblem to optimize MPS *)
If[verbose,
message2=PrintTemporary["Finding Q"];

];
time=AbsoluteTime[];
Q=FindQ[Heff,EnergyPerSite+shift];
	newEnergyPerSite="EnergyPerSite"/.Q;

time2=AbsoluteTime[];
findQtime=time2-time;
If[verbose,
NotebookDelete[message2];
message2=PrintTemporary["Finding Q took time: "<>ToString[findQtime]<>"s"];
];

(* Use SVD to decompose two-sites into 1*)
opts={"tolerance"->(tolerance//."tolerance"->0.),
"trunc"->(maxbond//."maxbond"->Min[Dimensions[Q]])
};
If[verbose,
message3=PrintTemporary["Splitting sites with SVD"];

];
time=AbsoluteTime[];
decomp=decompose[Q,Sequence@@opts];
	L="L"/.decomp;
	SingularValues="SingularValues"/.decomp;
time2=AbsoluteTime[];
SVDtime=time2-time;
If[verbose,
NotebookDelete[message3];
message3=PrintTemporary["Splitting sites with SVD took time: "<>ToString[SVDtime]<>"s"];
];
(* Construct new ZL and HL *)
time=AbsoluteTime[];
newZL=makeZL[L];
newHL=makeHL[L,HL,ZL,h,J,newEnergyPerSite];
HZLtime=AbsoluteTime[]-time;
(* Construct Expectation Values *)
Z=singlesiteexpect[L,{"tensor"->pauliZ},SingularValues];
X=singlesiteexpect[L,{"tensor"->pauliX},SingularValues];
ZZ=twositeexpect[L,{"tensor"->pauliZ},{"tensor"->pauliZ},SingularValues];
alltime=AbsoluteTime[]-time0;

If[verbose,NotebookDelete[{message1,message2,message3}]];

(* Return Result *)
Join[{"ZL"->newZL,"HL"->newHL,"E"->newEnergyPerSite,"J"->J,"h"->h,"tolerance"->tolerance,"maxbond"->maxbond,"L"->L,"singularvalues"->SingularValues,"iteration"->iteration+1,"Heff"->Heff,"Q"->Q,"guess"->EnergyPerSite+shift,"shift"->shift,"Z"->Z,"X"->X,"ZZ"->ZZ,"makeHtime"->makeHtime,"alltime"->alltime,"findQtime"->findQtime,"SVDtime"->SVDtime,"HZLtime"->HZLtime},If[verbose,{"verbose"->True},{}]]]

Verify that you get what I do (modulo the output of “HL” -- since I probably output more intermediate results than you want

i1=iterateIDMRG[{"ZL"->ZL1,"HL"->HL1,"J"->1.2,"h"->0.2,"E"->-10}]

i1=iterateIDMRG[{"ZL"->ZL1,"HL"->HL1,"J"->1.2,"h"->0.2,"E"->-10}]

Same thing but throw away the small singular values

i1s=iterateIDMRG[{"ZL"->ZL1,"HL"->HL1,"J"->1.2,"h"->0.2,"E"->-10,"tolerance"->10^-6}]

i1s=iterateIDMRG[{"ZL"->ZL1,"HL"->HL1,"J"->1.2,"h"->0.2,"E"->-10,"tolerance"->10^-6}]

Instead truncate to a fixed bond dimension

i1t=iterateIDMRG[{"ZL"->ZL1,"HL"->HL1,"J"->1.2,"h"->0.2,"E"->-10,"maxbond"->1}]

AbsoluteTiming[i2=iterateIDMRG[{"ZL"->ZL2,"HL"->HL2,"J"->1.2,"h"->0.2,"E"->-10,"verbose"->True}]]

Flow Control

Iterating by hand

The way we have set things up, we can iterate by hand -- which is useful for debugging

Here is an initial condition in what should be the ordered phase

start05={"ZL"->ZL2,"HL"->HL2,"E"->-10,"J"->1.,"h"->0.5,
"maxbond"->100,"shift"->-10.,"tolerance"->10^-6}

Here is running one iteration

run1=iterateIDMRG[naive05]

We have things set up to output a bunch of intermediate results -- which are there mostly for debugging.  In production code you would likely output less to save disk space / memory, maybe having a “verbose” argument that controls output.

We can get a list of all stored quantities by

run1[[All,1]]

And can look at them with commands like

MatrixForm[e1="tensor"/.("Heff"/.run1)]

We can run again with

run2=iterateIDMRG[run1];

"E"/.run2

Suppose we want to change the bond dimension.  We can do that with

run2b=iterateIDMRG[run1/.( "maxbond"->_)->("maxbond"->2)];

"E"/.run2b

"singularvalues"/.run2

"singularvalues"/.run2b

Built in flow control

Mathematica has some built in flow controls: Nest, NestList, FixedPoint, FixedPointList -- which allow one-liners that will do much of what we want.

For example, this will iterate at most 100 times, stopping when sequential energies differ by less than 10^-6, and sequential <Z> differ less than by 10^-8

AbsoluteTiming[run05a=FixedPointList[iterateIDMRG,start05,100,SameTest->((Abs[("E"/.#1)-("E"/.#2)]<10^-6)&&(Abs[("Z"/.#1)-("Z"/.#2)]<10^-8)&)];]

"E"/.run05a

"Z"/.run05a

"singularvalues"^2/.run05a[[-1]]

If we move closer to the critical point, the entanglement entropy increases

start09=start05 /.("h"->_)->("h"->0.9)

AbsoluteTiming[run09a=FixedPointList[iterateIDMRG,start09,100,SameTest->((Abs[("E"/.#1)-("E"/.#2)]<10^-6)&&(Abs[("Z"/.#1)-("Z"/.#2)]<10^-8)&)];]

ListPlot[{Legended["Z"/.run05a,"h=0.5 J"],Legended["Z"/.run09a,"h=0.9 J"]},AxesLabel->{"iteration","<Z>"},PlotRange->{0,1}]

ListPlot[entropy=-#.Log[#]&/@(("singularvalues"/.Rest[run09a])^2),AxesLabel->{"iteration","S"},PlotRange->All]

Custom Flow Control

While the simplicity of doing flow control with a single command like “FixedPointList” is convenient, it is often nice to write a custom iterator -- for example, which only saves the state at every jth iteration -- and also that updates the user about how long things are taking

Here is a typical custom loop

defaultparams={"maxiterations"->10,"accuracygoal"->10^-6,"zaccuracygoal"->10^-3,"reportinginterval"->1};

runIDMRG[DMRGparams_List,flowparams_List]:=
Module[
{maxiterations="maxiterations"/.flowparams//.defaultparams,
accuracygoal="accuracygoal"/.flowparams//.defaultparams,
zaccuracygoal="zaccuracygoal"/.flowparams//.defaultparams,
reportinginterval="reportinginterval"/.flowparams//.defaultparams,
h="h"/.DMRGparams,
abort=False,
result,i=0,run=True,DMRGstate=DMRGparams,newDMRGstate,
abortbutton,progressmessage,accuracy,lastE="E"/.DMRGparams,currentE="E"/.DMRGparams,zaccuracy,messageon=False,
buttonon=False,
time1,time2},
time1=AbsoluteTime[];
result=Reap[
While[run,
newDMRGstate=iterateIDMRG[DMRGstate];
accuracy=Abs[(lastE="E"/.DMRGstate)-(currentE="E"/.newDMRGstate)];
zaccuracy=Abs[("Z"/.DMRGstate)-("Z"/.newDMRGstate)];
DMRGstate=newDMRGstate;
i+=1;
If[Mod[i,reportinginterval]==0,
Sow[DMRGstate];
];
time2=AbsoluteTime[];
If[time2-time1>1 (*we have not written to the screen for over a second *),
time1=time2;
If[Not[buttonon],abortbutton=PrintTemporary[Button["Abort",abort=True]];buttonon=True];
If[messageon,
NotebookDelete[progressmessage]];
messageon=True;
progressmessage=
PrintTemporary["h="<>ToString[h]<>" E="<>ToString[currentE]<>" last E="<>ToString[lastE]<>
" i="<>ToString[i]<>" of "<>ToString[maxiterations]];
];
If[i>maxiterations,run=False];
If[(accuracy<accuracygoal)&&(zaccuracy<zaccuracygoal),run=False];
If[abort,run=False];
]];
NotebookDelete[abortbutton];
NotebookDelete[progressmessage];
Append[Flatten[result[[-1]],1],newDMRGstate]]

AbsoluteTiming[test1=runIDMRG[
{"ZL"->ZL2,"HL"->HL2,"E"->-10,"J"->1.,"h"->1.1,
"maxbond"->100,"shift"->-1.,"verbose"->False,"tolerance"->10^-6},
{"maxiterations"->100,"accuracygoal"->10^-6,"zaccuracygoal"->10^-3,"reportinginterval"->5}];]

ListPlot[{"iteration","E"}/.test1]

ListPlot[{"iteration","Z"}/.test1,PlotRange->All]

Looping over loops

Here is an example of running our iDMRG routines for a series of values of h, with the bond dimesion restricted to be 1.  This is almost our mean-field theory from lecture 1, except that the center two sites effectively have a slightly higher bond dimension

AbsoluteTiming[meanfielddat=
Module[{h=0,dh=0.1,state={"ZL"->ZL2,"HL"->HL2,"E"->-10,"J"->1.,"h"->0.,
"maxbond"->1,"shift"->-1.,"verbose"->False,"tolerance"->10^-6}},
Reap[While[h<3,
state=Last@runIDMRG[state,{"maxiterations"->100,"accuracygoal"->10^-6,"zaccuracygoal"->10^-3,"reportinginterval"->100,"progressinterval"->20}];
h=h+dh;
Sow[state];
state=state/.{("h"->_)->("h"->h),("iteration"->_)->("iteration"->1)}]][[-1,1]]];]

Dimensions[meanfielddat]

ListPlot[{Labeled[{"h","Z"}/.meanfielddat,"<Z>",0.5],Labeled[{"h","X"}/.meanfielddat,"<X>",0.5]},AxesLabel->{"h",None},PlotLabel->"Max bond dimension=1"]

Here is a function that goes through a list of states, and reruns the iDMRG solver with a new bond-dimension

refine[datset_,newmaxbond_]:=Module[{state=datset/.{("maxbond"->_)->("maxbond"->newmaxbond),("iteration"->_)->("iteration"->1)}},Last@runIDMRG[state,{"maxiterations"->100,"accuracygoal"->10^-6,"zaccuracygoal"->10^-3,"reportinginterval"->100}]]

AbsoluteTiming[lev2dat=refine[#,2]&/@meanfielddat;]

ListPlot[{{"h","Z"}/.meanfielddat,{"h","Z"}/.lev2dat}]

Near the phase boundaries it takes more iterations to converge.

ListPlot[{{"h","iteration"}/.meanfielddat,{"h","iteration"}/.lev2dat},PlotRange->All,AxesLabel->{"h","Number of Iterations"}]

Exercise 1

Run the iDMRG code with a max bond-dimension of 3, 5, 10.  Use these results to estimate the true location of the phase boundary.

Exercise 2

Sit at the phase boundary (h=1).  This is the place where the MPS ansatz is worst (as the correlation length diverges).  Starting from the state “start1” given below, iterate 100 times, plotting  <Z> vs iteration number.  Also plot correlation length vs iteration number, and entanglement entropy vs iteration number.  You should see that convergence  is extremely slow.

start1={"ZL"->ZL2,"HL"->HL2,"E"->-10,"J"->1.,"h"->1.,
"maxbond"->20,"shift"->-10.,"tolerance"->10^-6}