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

Jump to: navigation, search
(Spin-1/2 Chain)
(Spin-1/2 Chain)
Line 22: Line 22:
<math> \langle S^z_i S^z_j \rangle \sim (-1)^{|i-j|} \frac{\sqrt{\ln|i-j|}}{|i-j|} </math> .
<math> \langle S^z_i S^z_j \rangle \sim (-1)^{|i-j|} \frac{\sqrt{\ln|i-j|}}{|i-j|} </math> .
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 <math>-1/2</math>. There is also an additional square root of a logarithm 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.
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 <math>-1</math>. There is also an additional square root of a logarithm 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===
===Spin-1 Chain===

Revision as of 09:08, 17 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{\sqrt{\ln|i-j|}}{|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. There is also an additional square root of a logarithm 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


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

The most important correlation functions in many-body physics are two-point correlators, i.e. correlators that involve two sites i and j, such as \langle S^+_i S^-_j \rangle. Short-ranged ones determine energies (in the typical short-ranged Hamiltonians of correlation physics), long-ranged ones determine correlation lengths.

Another Go At The Energy Per Bond

As already mentioned above, the ground state energy per bond in both spin-1/2 and spin-1 chain is given by

e_0(i) = \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

This gives the energy of each bond individually, but we are interested in the thermodynamic limit, where all bonds are on equal footing and hence should have the same energy unless there is some physical breaking of translational invariance.

Obviously, the bonds that are closest to the thermodynamic limit behaviour are those in the chain center. So, the direct approach would be to calculate e_0(L/2) and extrapolate it first in D for fixed L and then in L.

Before you do this, however, plot for some values of D and not too small L e_0(i) versus i (as a check of the program, you may also consider the three contributions individually before you do the sum. What relationship between them should exist?).

What do you observe for spin-1? And what for spin-1/2?

For the spin-1/2 chain, bond energies oscillate strongly between odd and even numbered bonds. This is because the open ends make themselves felt very strongly due to criticality and because the spin-1/2 chain is on the verge of dimerization, i.e. a spontaneous breaking of translational symmetry of the ground state down to a periodicity of 2. It is therefore more meaningful to extrapolate the average energy of a strong and a weak bond; you immediately gain lots of accuracy. This is yet another example that it is worthwhile to have a close look at the actual output of DMRG by considering various local or (here) almost local observables.

Spin-Spin Correlations: Spin-1/2

Take a relatively long chain (say, L=192), and calculate  \langle S^z_i S^z_j \rangle for various increasing D.

Now plot C_l = \langle S^z_{L/2-l/2} S^z_{L/2+l/2} \rangle where you round the positions such that their distance is l. The purpose of this is to center the correlators about the chain center to make boundary effects as small as possible; there are also other ways of doing this (like averaging over several correlators with same site distance, also more or less centered). As we expect a power law, use a log-log plot. Take absolute values or multiply out the antiferromagnetic factor (-1)^l.

What you should see, is a power law on short distances, but a faster (in fact, exponential) decay for larger distances. This has two reasons: (i) the finite system size cuts off the power-law correlations; but as we took a large system size here, this should not matter too much. (ii) DMRG's algorithmic structure effectively generates correlators which are superpositions of up to D^2 purely exponential decays, and therefore can only mimic power laws by such superpositions - at large distance, the slowest exponential decay will survive all the others, replacing the power law by an exponential law. The larger you choose D, the further you push out this crossover.

Spin-Spin Correlations: Spin-1

In the spin-1 chain, we do expect exponential decay (with an analytic modification), so the exponential nature of the correlators of DMRG should fit well. Again, choose a long chain (say,L=192), and calculate  \langle S^z_i S^z_j \rangle for various increasing D.

Now plot C_l = \langle S^z_{L/2-l/2} S^z_{L/2+l/2} \rangle where you round the positions such that their distance is l, as before. As we expect an exponential law, use a log-lin plot, again eliminating the negative signs.

From the log-lin plot, extract a correlation length. It will depend (and in fact monotonically increase with) D. Has it converged when you reach e.g. D=300? How does this compare to the convergence for the same number of states of local or quasi-local observables such as magnetization or energy?

In fact, the calculation of correlation lengths is much harder to converge than that of the local quantities. This is due to the fact that a more profound algorithmic analysis reveals DMRG to be an algorithm geared especially well to the optimal representation of local quantities, not so much non-local ones as long-ranged correlators.

Sometimes There Is A Way Out

In the special case of the spin-1 chain, we have a loophole for the calculation of the correlation length, which is related to the weird observation that the first excitation was not a bulk excitation. It can be shown that a good toy model for a spin-1 chain is given as follows: at each site of a spin-1, you put two spin-1/2, and construct the spin-1 states from the triplet states of the two spin-1/2 at each site. The ground state is then approximated quite well by a state where you link two spin-1/2 on neighbouring sites by a singlet state.

In this construction, for open boundary conditions (but not periodic ones), on the first and on the last site there will be two lonely spins-1/2 without partner. These two spins-1/2 can form 4 states among themselves, which in the toy model are degenerate: the ground state is four-fold degenerate. In the real spin-1 chain, this four-fold degeneracy (from one state of total spin 0 and three of total spin 1) is only achieved in the thermodynamic limit when the two spins are totally removed from each other. This is why there was no gap between magnetization sectors 0 and 1. The first bulk excitation needs magnetization 2.

What we can do to cure this, is to attach one spin-1/2 before the first and after the last site, taking the same bond Hamiltonian, that link up to the two lonely spins by a singlet state. You may check that now the gap is between magnetization sectors 0 and 1!

In order to calculate the correlation length, one can also play the following trick: attach only one spin-1/2 at one end. This means that the ground state will now be doubly degenerate, in magnetization sectors +1/2 or -1/2, and be characterized by the boundary site where there is NO spin-1/2 attached carrying finite magnetization, that decays into the bulk, with the correlation length.

For a chain of length L=192 and D=200, calculate the ground state magnetization. Plot it (eliminating the sign oscillation) versus site in a log-lin plot and extract the correlation length. What do you get?