Difference between revisions of "ALPS 2 Tutorials:DMRG-01 DMRG"

From ALPS
Jump to: navigation, search
Line 77: Line 77:
 
ADRIAN, THIS SHOULD BE DONE BY YOU.
 
ADRIAN, THIS SHOULD BE DONE BY YOU.
  
Starting the code generally. How do we assert that enough sweeps have been done? What will be written here might affect a bit the later parts.
+
Starting the code generally. How do we assert that enough sweeps have been done? What will be written here might affect a bit the later parts. Good quantum number magnetization.
  
 
==Ground State Energies==
 
==Ground State Energies==
Line 126: Line 126:
 
Obviously, we have to be able to get access to the first excited state and its energy. DMRG fundamentally knows two ways of doing this, one which works always, but is not as neat, and another one, which is very clean, but does not work under all circumstances.
 
Obviously, we have to be able to get access to the first excited state and its energy. DMRG fundamentally knows two ways of doing this, one which works always, but is not as neat, and another one, which is very clean, but does not work under all circumstances.
  
# The pedestrian way
+
# The pedestrian way is to set up a DMRG calculation that calculates both states at the same time. However, for a given number of states the accuracy will somewhat decrease, as two different quantum states both have to be described well.
# The smarter way  
+
# The smarter way reduces the gap calculation to the calculation of two ground states. In many quantum systems, the ground state and the first excited state differ by a good quantum number and therefore are both ground states in the respective sectors. For example, for the spin-1/2 chain, the ground state is a singlet of total spin 0, and hence the ground state in the sector of magnetization 0. The first excited state is a triplet of total spin 1, i.e. consists of one excited state of magnetization 0, and the ground states of the sectors of magnetization +1 and -1 respectively. It can therefore be calculated as the ground state in magnetization sector +1.
 +
 
 +
Let us do this calculation first for the spin-1/2 chain. In a first attempt, fix <math>D=50,100,150</math> and calculate the gap for lengths <math>L=32,64,96,128</math>. For fixed <math>D</math>, plot the gap versus <math>1/L</math>. What you should see is that for small <math>D</math>, the results will not lie on a straight line passing through 0, but they will curve up from it. This behaviour gets better when <math>D</math> gets larger. Discuss why this might be.
 +
 
 +
In a second, more meaningful attempt, fix the lengths <math>L=32,64,96,128</math> and vary <math>D=50,100,150,200</math> in order to extrapolate the gap for each fixed length in <math>D</math> (or, as explained above, the truncation error). What does the plot of the gap versus <math>1/L</math> look like now?
  
 
===Extrapolating The Gap To The Thermodynamic Limit===
 
===Extrapolating The Gap To The Thermodynamic Limit===
  
What is the gap of the S=1 antiferromagnetic Heisenberg chain?
+
The case of the spin-1/2 chain is a bit frustrating, because all you will be able to say, even if you push the computer to its limits, is that the gap seems to be extremely small to the best of your abilities and therefore is likely to vanish. But who can tell you that you are not looking at a case where the gap is, say, <math>e^{-50}</math>? This of course is a sobering reminder of the limits of even a highly accurate numerical method.
 +
 
 +
Let us therefore turn to a more rewarding question, what is the gap of the spin-1 antiferromagnetic Heisenberg chain?
 +
 
 +
Here, there is a nasty twist, which we will at the moment only state and act on, but explain only later: Calculate the gap not between the ground states of the magnetization sectors 0 and 1, but 1 and 2. If you wish, do it also for 0 and 1, for later reference, but the following refers to 1 and 2.
  
Assume you have <math>\Delta (L)</math> with machine precision, either by a suitable extrapolation or by a very high accuracy calculation. If you don't want to do the former, calculate the gap for system sizes <math>L=8,16,32,48,64,96,128,192,256</math> with <math>D=300</math> states each and 5 sweeps.  
+
Assume you have <math>\Delta (L)</math> with machine precision, either by a suitable extrapolation as discussed above or by a very high accuracy calculation. If you don't want to do the former, calculate the gap for system sizes <math>L=8,16,32,48,64,96,128,192,256</math> with <math>D=300</math> states each and 5 sweeps.  
  
 
As the effects of the open ends will decrease as <math>1/L</math>, it always makes sense to
 
As the effects of the open ends will decrease as <math>1/L</math>, it always makes sense to
first plot the gaps <math>\Delta (L)</math> versus <math>1/L</math>. Produce such a plot.
+
first plot the gaps <math>\Delta (L)</math> versus <math>1/L</math>, as was already done in the spin-1/2 case. Produce such a plot.
  
What you see, is a curve that is quite straight for small L and then starts bending upward. What gap would you obtain if you extrapolate the  
+
What you see, is a curve that is quite straight for small L and then starts bending upward. What gap would you obtain if you extrapolate the linear part of the curve naively? (This question is relevant for situations where the correlation length of the chain is so long that it becomes hard to see the asymptotic behaviour on reachable length scales.) Is it over- or underestimated?
 +
 
 +
What gap do you read off if you take the longest chain you have? Is it over- or underestimated?
 +
 
 +
The ideal would be to have an idea of what the asymptotic behaviour (the curved part for long lengths) is like analytically to extrapolate. Do a plot of the gap as <math>\Delta (L)</math> versus <math>1/L^2</math>. What does the curve now look like for big lengths? Extrapolate the gap.
 +
 
 +
The last plot was in fact motivated by the following argument: from Haldane's analysis of the spin-1 chain by the nonlinear sigma model, one expects that the lowest lying excitations (which for periodic boundary conditions can be labeled by a momentum <math>k</math>) are around <math>k=\pi</math> and have an energy
 +
 
 +
<math>E(k) = E_0 + \sqrt{\Delta^2 + c^2 (k-\pi)^2}</math>.
 +
 
 +
For the open boundary conditions, we may approximate <math>k-\pi</math> by <math>1/L</math> (think about a particle in a box), which gives a finite-system size gap of
 +
 
 +
<math>\Delta(L) \approx \Delta \left( 1 + \frac{c^2}{2\Delta^2 L^2} \right) </math>
 +
 
 +
and indicates that in the asymptotic limit the convergence should essentially be as <math>1/L^2</math>. How close do you get to the result <math>\Delta=0.41052</math>?
 +
 
 +
For those, that also did the gap between the ground states of magnetisation sectors 0 and 1, show that the gap you get there is essentially zero. All others, take this result for granted and start worrying: why is the finite gap the right one and the vanishing gap the wrong one? Is this a physics lottery? In fact, there is a very good reason why the spin-1 chain shows this peculiar behaviour for open boundary conditions that can be found analytically; but even if we were not so fortunate as to know it, we could detect the problem right away! This can be done by the observation of local observables.
  
 
==Local Observables==
 
==Local Observables==
 +
 +
As local observables we consider observables that are linked to one specific site. In the case of spin chains, the meaningful local observable is the local magnetization <math>\langle S^z_i \rangle </math>.
 +
 +
===Excitations in the Spin-1 Chain===
 +
 +
Take a chain of length <math>L=96</math> and <math>D=200</math>. Calculate the local magnetization <math>\langle S^z_i \rangle </math> and plot it versus the site <math>i</math> for the ground states in the magnetisation sectors 0, 1, and 2.
 +
 +
What you should obtain is an essentially flat curve for sector 0, a magnetisation which is essentially concentrated at the chain ends for sector 1, and a magnetisation which is both at the chain ends and in the bulk of the chain for sector 2. This means that the first excitation of the open chain is a boundary excitation, which would not exist on a closed system, and the second excitation of the open chain is a boundary plus a bulk excitation, which is the one we are interested in. For an as of now unknown reason, the energy of the first bulk excitation therefore has to be extracted from comparing sectors 1 and 2.
 +
 +
The morale of the story is that by looking at this local observable, we can distinguish boundary from bulk excitations in the spin-1 chain.
 +
 +
===Magnetisation in the Spin-1/2 Chain ===
 +
 +
Repeat a similar calculation for the spin-1/2 chain in the lowest magnetisation sectors. What do you observe here?
  
 
==Correlation Functions==
 
==Correlation Functions==
 +
 +
  
 
==Sometimes There Is A Way Out ...==
 
==Sometimes There Is A Way Out ...==

Revision as of 19:51, 16 May 2010


Models: Heisenberg Spin Chains

For applications of DMRG, we consider two models, namely the spin-1/2 and the spin-1 antiferromagnetic Heisenberg chains of length L given by the following Hamiltonian,

H = \sum_{i=1}^{L-1} \frac{1}{2} (S^+_i S^-_{i+1} + S^-_i S^+_{i+1}) + S^z_i S^z_{i+1} .

The reason why we are choosing these two models, which you may already know from other tutorials, is that despite their superficial similarity they exhibit completely different physical behaviour and pose very different challenges to the DMRG algorithm. Let us briefly review their physical properties.

Spin-1/2 Chain

The ground state of the spin-1/2 chain can be constructed exactly by the Bethe ansatz; we therefore know its ground state energy exactly. In the thermodynamic limit L\rightarrow\infty the energy per site is given by

E_0 = 1/4 - \ln 2 = -0,4431471805599...

Ground state energies as such are of limited interest if not compared to other energies. But this one can serve as a beautiful benchmark of the DMRG method. Of more interest is whether the ground state is separated from the excited states by an energy difference that survives also in the thermodynamic limit, i.e. whether the gap is vanishing or not. For the spin-1/2 chain, the gap is 0.

At the same time, one may ask what the correlation between spins on different sites looks like. One knows for the infinitely long spin-1/2 chain that asymptotically (i.e. for |i-j| \rightarrow \infty)

 \langle S^z_i S^z_j \rangle \sim (-1)^{|i-j|} \frac{\ln|i-j|}{\sqrt{|i-j|}} .

This means that the spin-1/2 chain is critical, i.e. the antiferromagnetic correlations between spins decay with their distance following a power law; in this case the exponent of the power law is obviously -1/2. There is also an additional logarithmic correction which can be beautifully verified by DMRG calculations on very long chains, but given the very slow increase of the logarithm with its argument, we can ignore it in a first go.

Spin-1 Chain

For decades, people thought that the spin-1 chain would behave similarly, of course with some quantitative differences due to the different spin lengths. It came as a big surprise in 1982 when Duncan Haldane pointed out that there should be a fundamental difference between isotropic antiferromagnetic Heisenberg chains depending on the length of the spin, namely between half-integer spins (S=1/2,3/2,...) and integer spins (S=1), with the difference being most pronounced for small spin lengths. Hence, the spin-1 chain became the focus of strong interest, and in fact DMRG had some of its most important early applications right for this system.

Unlike the spin-1/2 chain, the spin-1 chain has no properties that can be calculated exactly by analytical means. We have to rely completely on numerics when it comes to quantitative statements.

The ground state energy per site is given by

 E_0 = -1.401484039 ... .

Again, the question of the existence of a gap is more important, and here one of the big differences to the spin-1/2 chain becomes visible: in the thermodynamic limit, the gap in the spin-1 chain is finite and given by

 \Delta = 0.41052

to five-digit accuracy.

The question for the behaviour of the spin-spin correlations leads to yet another big difference to the spin-1/2 case. The correlations read asymptotically (i.e. for |i-j| \rightarrow \infty)

 \langle S^z_i S^z_j \rangle \sim (-1)^{|i-j|} \frac{\exp (-|i-j|/\xi)}{\sqrt{|i-j|}} .

The dominant contribution is now the exponential decay which happens on a length scale \xi, the correlation length which in this particular case is found numerically to be \xi=6.02. There is an analytic (power law) correction by a square root of the distance in the denominator, but this is often neglected in calculations of the correlation length, as it is a slow contribution compared to the fast exponential decay. It would matter, of course, if the correlation length were much larger.

The spin-1 chain is therefore a prime example for a non-critical quantum system with finite gap and exponentially decaying correlations. As it will turn out, for DMRG this type of system is much easier to do.

Plan Of The Tutorial

What we want to achieve in the following tutorial, is to be able to calculate all the above quantities on our own using ALPS DMRG while learning about the principal pitfalls in this numerical project.

Vive la difference ...

The most important difference to other numerical methods is that DMRG prefers open boundary conditions, such that there are two chain ends at site 1 and L, not a closed loop as for example exact diagonalization and most analytical methods would prefer. This was already implicit in the notation of the Hamiltonian above and will lead to some of the more subtle aspects of DMRG calculations.

Running The Code

General Remarks

Before we start, let us briefly discuss the inner logic of the DMRG algorithm without discussing it in full detail. Given a one-dimensional quantum system with local state spaces of dimension d, where d=2S+1 for spins of length S, the Hilbert space dimension explodes exponentially as d^L with system size L. Exact diagonalization achieves exact results in this exponentially large Hilbert space, at the price of small system sizes. Quantum Monte Carlo gives approximate results by stochastically sampling this large space, reaching much larger system sizes. The density-matrix renormalization group (DMRG) tries yet another approach, namely to identify very small subspaces of the exponentially large Hilbert space which are hoped to contain good, very good, even excellent approximations to the states of interest such as the ground state.

A first key control parameter is therefore D, which controls the number of states in the subspace. DMRG is monotonic in this parameter: the larger it is, the larger is the subspace and the better the approximation can be. There is also an exact limit: if D\rightarrow d^L, no states are discarded and the solution would be exact. This is however of no practical relevance; if such a large number of states could be achieved on the computer, exact diagonalization is a superior alternative. A remark on notation: given DMRG history, D comes under various names, like matrix dimension or number of block states.

The second key control parameter is of course system size L.

The third control parameter(s) can only be understood by looking even closer at the DMRG algorithm. In order to find the best approximation to a state, DMRG proceeds in two steps:

  1. In a first step (so-called infinite-system DMRG) the algorithm tries to find good subspaces by iteratively analyzing chains of length 2, 4, 6, until the desired system size L is reached. The procedure consists of splitting the chain in every iteration and insert two new sites at the center; the name comes from the fact that this procedure can of course be carried on infinitely, to take L to infinity; but don't expect very meaningful results as you approach infinity! A second remark is that this procedure favours chains of even length for DMRG treatment.
  2. In a second step (so-called finite-system DMRG) DMRG deals with the fact that the subspace selection for shorter chains could not yet take into account all the quantum fluctuations and correlations that would be present in the chain of final length L. What the method does, is to go through a series of further iterations to improve the quality of the subspaces. One such iteration looking at all sites of a chain is referred to as a sweep in DMRG. The number of sweeps is the last important control parameter: if it is chosen to small, the precision reachable for a given D is not achieved; if it is chosen too large, calculational effort is wasted, although it is of course always good to err on the safe side.

In a last remark, let us consider the truncation error, which is a good indicator of the accuracy achieved by a DMRG run. In a simplified perspective, at each point in the algorithm DMRG makes one step in the direction of exponential growth of state space and then asks how much accuracy can be retained if not allowing that step, by means of an analysis of a density matrix regarding the distribution of weights (eigenvalues) corresponding to its eigenstates. The approximations of DMRG are then reflected in the fact that some statistical weight has to be discarded, which is the so-called truncation error. In many DMRG applications, it can be as small as 10^{-12}, showing that the approximations made by DMRG are extremely light, which is the reason for the enormous success of the method. For the purpose of the tutorial it is important to know that the error in local quantities (energies, magnetizations, ...) is roughly proportional to (but usually quite a bit larger than) the truncation error, provided we are converged in the number of sweeps.

The ALPS DMRG Code and Its Control Parameters

ADRIAN, THIS SHOULD BE DONE BY YOU.

Starting the code generally. How do we assert that enough sweeps have been done? What will be written here might affect a bit the later parts. Good quantum number magnetization.

Ground State Energies

The first question we usually ask is for the ground state | \psi_0 \rangle and its energy E_0. Here we have to distinguish two cases.

First, we might be interested in the ground state energy for a given Hamiltonian on a chain of a given length L. Secondly, we might be interested in the energy per site (or per bond) in the thermodynamic limit.

Fixed Length Ground State Energies

Consider chains of length L=32,64,96,128. Both for spin-1/2 and spin-1, set up ground state energy calculations for numbers of states D=50,100,150,200,300. For each length, tabulate the truncation error and the ground state energies as a function of D. Experiment carefully with the number of sweeps to assure that for a given length and number of states your result is actually converged.

  1. For each system length and spin length, try to establish the connection between the accuracy of the total energy and the truncation error by plotting total energy vs. truncation error.
  2. Observe how the convergence in D deteriorates with length for spin-1/2 and spin-1. Apart from a global factor of order of the length, do you see a difference between the convergence behaviour in the two cases? Hint: What you should see is, that but for the global factor, the convergence for large system sizes is only weakly dependent of length for the spin-1 chain, but much more strongly dependent for the spin-1/2 chain. This is because the spin-1 chain physics is dominated by segments of length of the correlation length, whereas for the spin-1/2 chain there is no finite length scale because of criticality.
  3. Try to extrapolate ground state energies for each chain length to the D\rightarrow\infty limit.

Ground State Energies Per Site (Bond)

If we look closely at the Hamiltonian, the energy of a chain of length L does not sit on the L sites, but on the L-1 bonds. A first (naive) attempt therefore consists in taking the results of the last simulations and calculating

e_0 = \frac{E_0(L)}{L-1}.

Do you get really close to the values listed in the introduction? What is the wrong underlying assumption?

The correct way is to eliminate the effect of the open boundary conditions by considering the energy of one bond at the center of the chain. There are two ways of doing it.

  1. Calculate the ground state energy of two chains of length L and L+2, again for the lengths already mentioned above, and calculate e_0 = (E_0(L+2) - E_0 (L))/2 as the energy per bond. What do the results look like now?
  2. The less costly and usual way would be to use correlators (as discussed further below, so postpone this exercise until then) between neighbouring sites and use

e_0 = \frac{1}{2} \langle S^+_i S^-_{i+1}\rangle  + \langle S^-_i S^+_{i+1}\rangle ) + \langle S^z_i S^z_{i+1} \rangle

for sites i,i+1 at the chain center.

Calculating The Gap

As already mentioned, the energy gap of a quantum system is given by the energy difference between the first excited state and the ground state

\Delta = E_1 - E_0

in the thermodynamic limit. This means we have to solve two problems, (i) the calculation of

\Delta(L) = E_1 (L) - E_0 (L)

for finite system sizes and (ii) extrapolate \Delta (L) to the thermodynamic limit L= \infty. The latter is not specific to DMRG, but because of its preference for open boundary conditions somewhat more complicated than in the more usual case of periodic boundary conditions.

Getting The Gap For Finite Systems

Obviously, we have to be able to get access to the first excited state and its energy. DMRG fundamentally knows two ways of doing this, one which works always, but is not as neat, and another one, which is very clean, but does not work under all circumstances.

  1. The pedestrian way is to set up a DMRG calculation that calculates both states at the same time. However, for a given number of states the accuracy will somewhat decrease, as two different quantum states both have to be described well.
  2. The smarter way reduces the gap calculation to the calculation of two ground states. In many quantum systems, the ground state and the first excited state differ by a good quantum number and therefore are both ground states in the respective sectors. For example, for the spin-1/2 chain, the ground state is a singlet of total spin 0, and hence the ground state in the sector of magnetization 0. The first excited state is a triplet of total spin 1, i.e. consists of one excited state of magnetization 0, and the ground states of the sectors of magnetization +1 and -1 respectively. It can therefore be calculated as the ground state in magnetization sector +1.

Let us do this calculation first for the spin-1/2 chain. In a first attempt, fix D=50,100,150 and calculate the gap for lengths L=32,64,96,128. For fixed D, plot the gap versus 1/L. What you should see is that for small D, the results will not lie on a straight line passing through 0, but they will curve up from it. This behaviour gets better when D gets larger. Discuss why this might be.

In a second, more meaningful attempt, fix the lengths L=32,64,96,128 and vary D=50,100,150,200 in order to extrapolate the gap for each fixed length in D (or, as explained above, the truncation error). What does the plot of the gap versus 1/L look like now?

Extrapolating The Gap To The Thermodynamic Limit

The case of the spin-1/2 chain is a bit frustrating, because all you will be able to say, even if you push the computer to its limits, is that the gap seems to be extremely small to the best of your abilities and therefore is likely to vanish. But who can tell you that you are not looking at a case where the gap is, say, e^{-50}? This of course is a sobering reminder of the limits of even a highly accurate numerical method.

Let us therefore turn to a more rewarding question, what is the gap of the spin-1 antiferromagnetic Heisenberg chain?

Here, there is a nasty twist, which we will at the moment only state and act on, but explain only later: Calculate the gap not between the ground states of the magnetization sectors 0 and 1, but 1 and 2. If you wish, do it also for 0 and 1, for later reference, but the following refers to 1 and 2.

Assume you have \Delta (L) with machine precision, either by a suitable extrapolation as discussed above or by a very high accuracy calculation. If you don't want to do the former, calculate the gap for system sizes L=8,16,32,48,64,96,128,192,256 with D=300 states each and 5 sweeps.

As the effects of the open ends will decrease as 1/L, it always makes sense to first plot the gaps \Delta (L) versus 1/L, as was already done in the spin-1/2 case. Produce such a plot.

What you see, is a curve that is quite straight for small L and then starts bending upward. What gap would you obtain if you extrapolate the linear part of the curve naively? (This question is relevant for situations where the correlation length of the chain is so long that it becomes hard to see the asymptotic behaviour on reachable length scales.) Is it over- or underestimated?

What gap do you read off if you take the longest chain you have? Is it over- or underestimated?

The ideal would be to have an idea of what the asymptotic behaviour (the curved part for long lengths) is like analytically to extrapolate. Do a plot of the gap as \Delta (L) versus 1/L^2. What does the curve now look like for big lengths? Extrapolate the gap.

The last plot was in fact motivated by the following argument: from Haldane's analysis of the spin-1 chain by the nonlinear sigma model, one expects that the lowest lying excitations (which for periodic boundary conditions can be labeled by a momentum k) are around k=\pi and have an energy

E(k) = E_0 + \sqrt{\Delta^2 + c^2 (k-\pi)^2}.

For the open boundary conditions, we may approximate k-\pi by 1/L (think about a particle in a box), which gives a finite-system size gap of

\Delta(L) \approx \Delta \left( 1 + \frac{c^2}{2\Delta^2 L^2} \right)

and indicates that in the asymptotic limit the convergence should essentially be as 1/L^2. How close do you get to the result \Delta=0.41052?

For those, that also did the gap between the ground states of magnetisation sectors 0 and 1, show that the gap you get there is essentially zero. All others, take this result for granted and start worrying: why is the finite gap the right one and the vanishing gap the wrong one? Is this a physics lottery? In fact, there is a very good reason why the spin-1 chain shows this peculiar behaviour for open boundary conditions that can be found analytically; but even if we were not so fortunate as to know it, we could detect the problem right away! This can be done by the observation of local observables.

Local Observables

As local observables we consider observables that are linked to one specific site. In the case of spin chains, the meaningful local observable is the local magnetization \langle S^z_i \rangle .

Excitations in the Spin-1 Chain

Take a chain of length L=96 and D=200. Calculate the local magnetization \langle S^z_i \rangle and plot it versus the site i for the ground states in the magnetisation sectors 0, 1, and 2.

What you should obtain is an essentially flat curve for sector 0, a magnetisation which is essentially concentrated at the chain ends for sector 1, and a magnetisation which is both at the chain ends and in the bulk of the chain for sector 2. This means that the first excitation of the open chain is a boundary excitation, which would not exist on a closed system, and the second excitation of the open chain is a boundary plus a bulk excitation, which is the one we are interested in. For an as of now unknown reason, the energy of the first bulk excitation therefore has to be extracted from comparing sectors 1 and 2.

The morale of the story is that by looking at this local observable, we can distinguish boundary from bulk excitations in the spin-1 chain.

Magnetisation in the Spin-1/2 Chain

Repeat a similar calculation for the spin-1/2 chain in the lowest magnetisation sectors. What do you observe here?

Correlation Functions

Sometimes There Is A Way Out ...