• View
• Edit
• Rename
• History
• Print
• Backlinks

IDMRGIntroduction

In this tutorial, we will use infinite-size DMRG to obtain the ground-state energy of the spin-1/2 Heisenberg model, and some properties such as the correlation length.

Firstly, we need to obtain a lattice file for our system. There are predefined models for spin chains with {$SU(2)$}, {$U(1)$}, {$Z_2$} or no symmetry, which are named spinchain-su2, spinchain-u1, spinchain-z2 and spinchain, respectively. The differences are in the allowed interactions, since the Hamiltonian must respect the chosen symmetry. So spinchain-u1 requires that the Hamiltonian (and observables) conserve the z-component of spin, while spinchain-su2 requires that the Hamiltonian and observables preserve the full spin vector {$\vec{s}$}. In this tutorial, we will start without using any symmetry, using the spinchain model, however all of the variants contain the basic operators for isotropic spin-spin interactions, so spinchain-u1 or spinchain-su2 could be used instead. Indeed, this is a useful exercise, to see the effect of different symmetries.

To construct the lattice file:

spinchain -o lattice

We have used the default value of the spin (0.5), otherwise we could set it explicitly, for example a spin-3/2 chain we can add the option -S 1.5.

To obtain the groundstate wavefunction:

We use the mp-idmrg-s3e program. The Hamiltonian we want is lattice:H_J1, which represents nearest-neighbor isotropic spin exchange. This will be the option -H lattice:H_J1 to the idmrg program.

The first time we run the program, we need to specify that we want to create a new wavefunction (--create), and also specify the unit cell to be two sites (-u 2). We need this because the groundstate of the spin-1/2 chain has momentum {$\pi$}, so it is translationally invariant under a 2-site shift. If we were using {$U(1)$} or {$SU(2)$} symmetry we would also need to specify the quantum number per unit cell (option (-q 0), for a singlet state), but we don't need this for the model with no symmetries.

We also need to specify the number of states, and the number of sweeps. For our first calculation we are going to use m=50 states, and 200 sweeps should be enough for convergence.

mp-idmrg-s3e -w psi -H lattice:H_J1 -m 50x200 --create -u 2

This should only take a few seconds to run. At the end of the calculation, the wavefunction is explicitly orthogonalized to be translationally invariant, before it is saved to disk.

For any wavefunction, we can obtain some general information about it using the mp-info command. This displays various information, including the type of wavefunction (eg finite, infinite, or some more exotic options), symmetry and quantum number, the unit cell size, and more. It also shows the last entry of the history log, of commands that have been used that modified the wavefunction. (If you want the complete history log, use the mp-history command.) mp-info can also show much more information, including the bipartite entropy at each bond, the dimension of the MPS at each site, the density matrix and entangement spectrum, etc. Run mp-info with no parameters to see a help message and link to the online documentation.

In order to verify that the calculation has converged, we could run this command again. Since we want to start from the existing groundstate, we must remember to remove the --create option, otherwise it would overwrite the wavefunction with a new one (in fact the toolkit will refuse to do that, unless you also supply the --force option, to force overwritiing the file). Since we are starting from an existing wavefuntion, we don't need to supply the unit cell option -u either, and we also don't need to supply the Hamiltonian -- there is a default Hamiltonian stored as a wavefunction attribute, which you can see using the command mp-attr psi). So you can continue the DMRG calculation using

mp-idmrg-s3e -w psi -m 50x200

The resulting energy is the energy per unit cell. Since we used a unit cell of 2 sites, this is double the energy per site. The exact energy per site in the thermodynamic limit should be -ln(2) + 0.25. For the 2-site unit cell and no symmetry this this is approximately -0.88628383... This isn't the optimal energy for m=50 states because it uses a single-site algorithm with density matrix mixing, and the mixing perturbs the density matrix, and hence the basis of the MPS. If you want to get closer to the optimal state for a given basis size, see the 'Advanced Techniques' section below.

NOTE: The energy that is displayed on the screen while mp-idmrg-s3e is running isn't a true variational energy -- it represents a local energy of a wavefunction that isn't perfectly translationally invariant. So it that might be lower than the groundstate energy, but typically it is a bit higher. To get the true variational energy use the mp-imoments command:

mp-imoments psi

You can use this tool to calculate the expectation value of any triangular MPO -- by default it will evalute the expectation value of the Hamiltonian attribute of the wavefunction, but you can calculate some other expectation value by using mp-imoments <psi> <operator>.

The output of mp-imoments is in the form of a polynomial. The #moment column is the total degree of the polynomial, and each row gives the coefficient of the term at the given #degree. The total energy is a 1st degree polynomial, since it is a linear function of the system size. So we can interpret the output as {$E = -0.8862838 L$}, where {$L$} is the number of unit cells in the system (in the large {$L$} limit), or equivalently, as the energy per unit cell as {$E/L = -0.88628383$}.

It is often convenient to scale the energy in a different way to the wavefunction unit cell. Since the Hamiltonian is translationally invariant under a single site shift, typically energies of the Heisenberg model are reported as an energy per site, rather than per 2-site unit cell. If we want to display the output of mp-imoments scaled to a different size unit cell, we can supply the -u (or --unitcell) parameter. For example,

mp-imoments -u 1 psi

displays the energy per site.

We can also calculate the variance per unit cell, from the expectation value of {$\langle H^2 \rangle$}. To do this, we add --power 2, to calculate the second power of the operator. (We could equivalently use the operator lattice:H_J1^2).

mp-imoments psi --power 2

The expectation value {$\langle H^2 \rangle = E_0^2 L^2 + \sigma^2 L$}, where {$E_0$} is the energy per unit cell, and {$\sigma^2$} is the variance per unit cell. So the variance is given by the degree 1 part of the 2nd moment, which should be around 1E-5 in this example.

Observables

The main tool for calculating observables is mp-expectation. This takes an operator that acts on a finite section of lattice (which can be arbitrarily large, but cannot have infinite support), and determines the expectation value. For example, the nearest-neighbor correlation

mp-expectation psi lattice:"Sz(0) * Sz(1)"

For obtaining correlation functions, there is a short-cut tool mp-icorrelation. This tool is much more efficient for calculating long correlation functions than running mp-expectation lots of times.

For infinite systems, transfer matrix methods can be much more powerful than looking at real-space correlations. The basic tool for examining the transfer matrix is mp-ispectrum. This tool has lots of uses but its basic function is to calculate the spectrum of possible correlation lengths of a wavefunction. use it like

mp-ispectrum psi

By default, mp-ispectrum displays the 10 largest eigenvalues in each symmetry sector, and the corresponding correlation length. The transfer matrix is block-diagonal with respect to the global symmetries of the wavefunction. This means that the eigenspectrum are labelled by the quantum numbers of the system. Since we are not using symmetries yet, it shows all quantum numbers as 0.

Note that for an MPS, the correlation length is always finite. This is a property of the MPS ansatz itself. However, we know from the Bethe Ansatz that the spin-1/2 chain is quantum critical, with gapless excitations corresponding to a {$c=1$} Conformal Field Theory and power-law correlations corresponding to an infinite correlation length. How then can we reconcile this with a finite correlation length in our MPS calculation? The answer is that the correlation length scales with the basis size of the MPS. If you construct a plot of correlation length {$\xi$} versus basis size {$m$} then you can see this explicitly. For details on this, see the HOWTO at HOWTO.CorrelationLengthScaling?.

Using symmetries

If the Hamiltonian has global symmetries such as particle number conservation, or spin rotational invariance, it is generally advantageous to use these symmetries in the calculation. Preserving a symmetry enforces a block structure into the MPS, which gives a useful speedup, which in some cases can be many orders of magnitude. It also enforces that the obtained groundstates exactly preserve the symmetry. This is more important for iMPS than for finite MPS, since variationally an iMPS will tend to break global symmetries in order to obtain a lower variational energy, even if they are symmetries that would not be broken in the exact groundstate.

We can demonstrate that the wavefunction we obtained above breaks various symmetries by calculating observables and other properties.

It breaks translational invariance and spin rotation symmetry because the expectation value of {$Sz$} is not zero. eg,

mp-expectation psi lattice:"Sz(0)"
mp-expectation psi lattice:"Sz(1)"
mp-expectation psi lattice:"Sp(0)"

The expectation value of {$S^+$} should always be zero for a state that preserves the z-component of spin, but here we see that the value is relatively big. The total Sz operator {$\sum_i S^z(i)$} is also non-zero, although it is small it is much bigger than a rounding error. Since the wavefunction is translationally invariant under a 2-site shift you can calculate the total Sz per unit cell as mp-expectation psi lattice:"Sz(0) + Sz(1)", but an alternative way of calculating this operator is using mp-imoments:

mp-imoments psi lattice:"sum_unit(Sz(0))"

The sum_unit function constructs a triangular MPO, as the sum of all translations of the lattice unit cell of the given operator. So sum_unit(Sz(0)) is the operator ({$\sum_i S^z(i)$}.

{$U(1)$} symmetry

We can construct a wavefunction that preserves the z-component of spin by using the spinchain-u1 model.

spinchain-u1 -o lattice-u1
mp-idmrg-s3e -w psi-u1 -H lattice-u1:H_J1 -m 50x200 --create -u 2 -q 0

We used a different name for the lattice and wavefunction, to avoid overwriting the files we constructed previously. We used a lattice constructed from spinchain-u1, and it is now necessary to specify the -q 0 option, meaning that the quantum number per unit cell is 0. (If you wanted to construct a magnetized excited state, you could set the quantum number per unit cell to be something different; -u 2 -q 1 would be a fully polarized ferromagnet; -u 4 -q 1 would be a 50\% polarised ferromagnet etc).

You will notice that this command runs even faster than the previous calculation that didn't use symmetry. This is because the MPS matrices now have a block structure, so the cost of operations on the matrices is reduced. This effect is more noticeable on larger scale calculations with more states kept.

If we calculate the total Sz operator for this state, we should see that it is now numerically zero, meaning it is smaller than about {$10^{-16}$}, which is smaller than the number of digits available in a double-precision floating point variable. Because we are using U(1) symmetry, the expectation value of a raising or lowering operator such as Sp(0) is identically zero (depending on the version of the toolkit you are using, attempting to calculate an expectation value of an operator that is not a rotational invariant will either give 0, or quit with an error because the quantum numbers to not match). We cannot calculate expectation values involving Sx or Sy using the U(1) symmetric model, because those operators do not transform irreducibly under U(1) rotations. That is, all operators must have a commutation relation with the Sz operator of

{$ [Sz, O] = nO$}

where {$n$} is a (half) integer, which is the label of the group representation. So {$n=0$} for the Sz operator, {$n=+1$} for the Sp operator, {$n=-1$} for the Sm operator, and so on. The Sx and Sy operators are not irreducible, so any expectation value using these operators needs to be re-expressed in terms of Sp, Sm, Sz before we can evaluate it.

If you use mp-ispectrum to look at the transfer matrix eigenvalues, you will notice that these are now split into many different quantum number sectors. This quantum number corresponds to the nature of the operator of the corresponding correlation function. That is, in sector [q], the correlation length corresponds to an operator that has that same quantum number. For example, the q=0 sector is scalar operators, something like a Sz-Sz or dimer-dimer correlation. The q=1 sector is the raising operator, so correlation lengths in this sector correspond to correlation functions of the form {$\langle S^+(0) S^-(x) \rangle$}, and so on.

{$SU(2)$} symmetry

We can construct a wavefunction for the spin-1/2 Heisenberg model that explicitly preserves {$SU(2)$} symmetry using the commands

spinchain-su2 -o lattice-su2
mp-idmrg-s3e -w psi-su2 -H lattice-su2:H_J1 -m 50x200 --create -u 2 -q 0

You should find that this takes approximately the same amount of time to run as it did for the {$U(1)$} symmetric code. Why is this? The group {$SU(2)$} is much bigger than {$U(1)$}, so we might expect that the efficiency should increase a lot, corresponding to the larger symmetry group. Well, if we compare the resulting energy, we notice that using {$SU(2)$} it is quite a bit better. Your exact numbers will vary, but typically with m=50 the {$SU(2)$} calculation will give an energy per site of around -0.4431468, compared with -0.4431419 using {$U(1)$}. The reason for this is that {$SU(2)$} is a non-abelian symmetry, meaning that the representations of the symmetry group are, in general, matrices of dimension > 1, rather than pure numbers. This corresponds to elements of the Hilbert space being multiplets, rather than single vectors. For a state in {$SU(2)$} with quantum number j, the dimension of the multiplet is {$2j+1$}. The effect of this is easiest to see by looking at the density matrix. We can see the density matrix at the edge of the unit cell (partition 0) using the command

mp-info psi-su2 -d -p 0

The important column is the 3rd column, #Degen, which is the degeneracy of the state (another name for the dimension of the representation). You will see that this corresponds to 2j+1 for quantum number sector j. This result is that the effective basis size is much greater than 50! Your results may vary slightly, but you should see that the total degree printed at the bottom (corresponding to the total effective number of states) is around 3 times bigger than the m=50 states that we requested. So although we only have m=50 states in total, the effective number of states is around 150! This is why the energy has improved compared with the {$U(1)$} code.

Calculating expectation values using {$SU(2)$} symmetry is more difficult, because we need to construct rotational invariant forms, and the algebra of {$SU(2)$} invariant operators is unfamiliar. The only local operator we have (apart from the identity) is the spin vector operator S itself, so we need to construct rotational invariant operators from combinations of spin operators. Something to keep in mind when manipulating spin vector operators is that they are axial vectors, meaning they are anti-hermitian: {$S^\dagger = -S$}. The toolkit defines some operations inner, dot and cross, which make it easier to construct rotational invariant operators. The difference between inner and dot for vector operators is inner(A,B) is equivalent to {$\vec{A}^\dagger \cdot \vec{B}$}, whereas dot(A,B) is equivalent to {$\vec{A} \cdot \vec{B}$}. Since the spin operator is anti-hermitian, in practice it is easier to use inner. For example,

mp-expectation psi-su2 lattice-su2:"inner(S(0), S(1))"

is the nearest-neighbor spin correlation function.

Advanced techniques

If we really want to push to get a well-converged wavefunction, we would do this several more times, and also decrease the density matrix mixing factor each time. The --mix-factor is the parameter used in the single-site subspace expansion algorithm http://arxiv.org/abs/1501.05504. This defaults to 0.02, which is good for early stages of the calculation, but for the best possible wavefunction at fixed basis size we want to slowly decrease this. For example,

mp-idmrg-s3e -w psi -H lattice:H_J1 -m 50x200 --mix-factor 1e-5
mp-idmrg-s3e -w psi -H lattice:H_J1 -m 50x200 --mix-factor 1e-7
mp-idmrg-s3e -w psi -H lattice:H_J1 -m 50x200 --mix-factor 0