Recent Changes - Search:

HomePage

PmWiki

pmwiki.org

Main /

MatlabDMRG

There is a simple Matlab program, downloadable in source form here.

There is a Python version of this code, described at Tutorials.SimpleDMRG. The Python page has a more extensive description of the code.

This is a simple program that reproduces the original implementation of DMRG. It is missing many features that an efficient DMRG code would contain.

  • The Superblock Hamiltonian doesn't need to be explicitly constructed. The Hamiltonian diagonalization in this code is {$O(m^4)$}, however it can be done in {$O(m^3)$} with a bit more care.
  • There is no initial guess vector for the Hamiltonian diagonalization. This means that at each iteration it is necessary to start from a random vector, even when the calculation has converged
  • The algorithm should converge to an infinite-size translationally-invariant MPS, but it is not obvious how to obtain it
  • There is no facility to calculate observables
  1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  2. %               DMRG for the  1D Heisenberg Model
  3. %               After S.R.White, Phys. Rev. Lett. 69, 2863 (1992)
  4. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  5. clear all
  6. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  7. %               Parameters
  8. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  9.  
  10. % Number of states kept m
  11. m=10;
  12. % Number of iterations.  Final lattice size is 2*Niter + 2
  13. NIter = 100;
  14. % exact energy, for comparison
  15. ExactEnergy = -log(2) + 0.25;
  16.  
  17. fprintf('Iter\tEnergy\t\tBondEnergy\tEnergyError\tTrunc\n');
  18. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  19. %               Intialize local operators
  20. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  21.  
  22. I= eye(2);
  23. Sz = [1/2 0;0 -1/2];
  24. Sp = [0 0;1 0];
  25. Sm = [0 1 ;0 0];
  26.  
  27. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  28. %               Initial blocks
  29. %               We assume reflection symmetry so we only need 1 block
  30. %               The operator acts on the inner-most site of the block
  31. %               +---------+    +---------+
  32. %               |        *|    |*        |
  33. %               +---------+    +---------+
  34. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35.  
  36. BlockSz = Sz;
  37. BlockSp = Sp;
  38. BlockSm = Sm;
  39. BlockI  = I;
  40. BlockH  = zeros(2);
  41. Energy = -0.75;   % energy of 2-site lattice
  42.  
  43. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  44. %               Begin main iterations
  45. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  46.  
  47. for l = 1:NIter
  48.  
  49.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  50.     %           Get the 2m-dimensional operators for the block + site
  51.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  52.  
  53.     BlockH = kron(BlockH, I) + kron(BlockSz, Sz) ...
  54.         + 0.5 * ( kron(BlockSp, Sm) + kron(BlockSm, Sp) );
  55.     BlockSz = kron(BlockI, Sz);
  56.     BlockSp = kron(BlockI, Sp);
  57.     BlockSm = kron(BlockI, Sm);
  58.     BlockI  = kron(BlockI, I);
  59.  
  60.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  61.     %                   HAMILTONIAN MATRIX forsuperblock
  62.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  63.  
  64.     H_super = kron(BlockH, BlockI) + kron(BlockI, BlockH) ...
  65.         + kron(BlockSz, BlockSz) ...
  66.         + 0.5 * ( kron(BlockSp, BlockSm) + kron(BlockSm, BlockSp) );
  67.  
  68.     H_super = 0.5 * (H_super + H_super');  % ensure H is symmetric
  69.  
  70.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  71.     %                   Diagonalizing the Hamiltonian
  72.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  73.  
  74.     LastEnergy = Energy;
  75.     opts.disp = 0;   % disable diagnostic information in eigs
  76.     opts.issym = 1;
  77.     opts.real = 1;
  78.     [Psi Energy] = eigs(H_super,1,'SA', opts);
  79.     EnergyPerBond = (Energy - LastEnergy) / 2;
  80.  
  81.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  82.     %                   Form the reduced density matrix
  83.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  84.  
  85.     [nr,nc]=size(Psi);
  86.     Dim = sqrt(nr);
  87.     PsiMatrix = reshape(Psi, Dim, Dim);
  88.     Rho = PsiMatrix * PsiMatrix';
  89.  
  90.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  91.     %                   Diagonalize the density matrix
  92.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  93.  
  94.     [V D] = eig(Rho);
  95.     [D Index] = sort(diag(D), 'descend');  % sort eigenvalues descending
  96.     V = V(:,Index);                     % sort eigenvectors the same way
  97.  
  98.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  99.     %                   Construct the truncation operator
  100.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  101.  
  102.     NKeep = min(size(D, 1), m);
  103.     T = V(:, 1:NKeep);
  104.     TruncationError = 1 - sum(D(1:NKeep));
  105.  
  106.     fprintf('%d\t%f\t%f\t%f\t%f\n',l, Energy, EnergyPerBond, ...
  107.         ExactEnergy-EnergyPerBond, TruncationError);
  108.  
  109.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  110.     %  Transform the block operators into the truncated basis
  111.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  112.  
  113.     BlockH  = T'*BlockH*T;
  114.     BlockSz = T'*BlockSz*T;
  115.     BlockSp = T'*BlockSp*T;
  116.     BlockSm = T'*BlockSm*T;
  117.     BlockI = T'*BlockI*T;
  118.  
  119. end
Get Code
Edit - History - Print - Recent Changes - Search
Page last modified on December 14, 2017, at 09:04 PM