Recent Changes - Search:

HomePage

PmWiki

pmwiki.org

Tutorials /

SimpleDMRG

(redirected from APCTPWorkshop.SimpleDMRG)

This is a simple DMRG code in Python. It follows quite closely the original implementation of the infinite system variant of DMRG from Steven R. White, Phys. Rev. Lett. 69, 2863 (1992). This is a Python port of a code first written in Matlab. It should run with python2 or python3. It requires the NumPy and SciPy libraries. You can get the Python source here

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. (Finding such an initial guess vector is the main development behind the PWFRG and iDMRG algorithms)
  • 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

This code makes use of reflection symmetry so that we only need to construct one set of block operators. If we didn't make use of this reflection symmetry, it would be tempting to add a site to left block as kron(block,site), and add a site the right block as kron(site,block), to give the usual DMRG form (block,site,site,block). In fact, it makes no difference whether we add sites to the block using kron(block,site) or kron(site,block), they differ simply by a unitary basis transformation which has no effect on the density matrix. Note that we could make more use of reflection symmetry, by requiring that the wavefunction is actually an eigenstate of reflection. This corresponds to requiring that the wavefunction in the Schmit basis (named PsiMatrix in the code) is a symmetric matrix. The code doesn't make use of this, but as an exercise, you can verify that PsiMatrix is indeed numerically symmetric.

  1. #####################################################
  2. # Simple DMRG example for the 1D Heisenberg model   #
  3. # Ian McCulloch, August 2017                        #
  4. # After S.R.White, Phys. Rev. Lett. 69, 2863 (1992) #
  5. #####################################################
  6.  
  7. import numpy as np
  8. import scipy
  9. import scipy.sparse.linalg
  10. import math
  11.  
  12. ##
  13. ## Initial parameters
  14. ##
  15.  
  16. # number of states kept m
  17. m = 10
  18. # number of iterations.  Final lattice size is 2*NIter + 2
  19. NIter = 100
  20.  
  21. # Exact energy per site, for comparison
  22. ExactEnergy = -math.log(2) + 0.25
  23.  
  24. print("  Iter  Size     Energy        BondEnergy   EnergyError   Truncation")
  25.  
  26. ##
  27. ## local operators
  28. ##
  29.  
  30. I = np.mat(np.identity(2))
  31. Sz = np.mat([[0.5,  0  ],
  32.              [0  , -0.5]])
  33. Sp = np.mat([[0, 0],
  34.                [1, 0]])
  35. Sm = np.mat([[0, 1],
  36.              [0, 0]])
  37.  
  38. ##
  39. ## Initial block operators.  At the start
  40. ## these represent a single site.  We assume
  41. ## reflection symmetry so we only need one block
  42. ## we can use on both the left and the right
  43. ##
  44.  
  45. BlockSz = Sz
  46. BlockSp = Sp
  47. BlockSm = Sm
  48. BlockI = I
  49. BlockH = np.zeros((2,2))  # Hamiltonian for 1-site system
  50.  
  51. Energy = -0.75        # initial energy for 2 sites
  52.                       # (we start iterations from 4 sites)
  53.  
  54. ##
  55. ## Begin main iterations
  56. ##
  57.  
  58. for i in range(0,NIter):
  59.  
  60.     ## Add a site to the block
  61.     BlockH = np.kron(BlockH, I) + np.kron(BlockSz, Sz) + \
  62.              0.5 * (np.kron(BlockSp, Sm) + np.kron(BlockSm, Sp))
  63.     BlockSz = np.kron(BlockI, Sz)
  64.     BlockSp = np.kron(BlockI, Sp)
  65.     BlockSm = np.kron(BlockI, Sm)
  66.     BlockI = np.kron(BlockI, I)
  67.  
  68.     ## 'Superblock' Hamiltonian
  69.  
  70.     H_super = np.kron(BlockH, BlockI) + np.kron(BlockI, BlockH) + \
  71.               np.kron(BlockSz, BlockSz) + 0.5 * (np.kron(BlockSp, BlockSm) + \
  72.                                                  np.kron(BlockSm, BlockSp))
  73.  
  74.     ## Diagonalize the Hamiltonian
  75.  
  76.     LastEnergy = Energy
  77.     E, Psi = scipy.sparse.linalg.eigsh(H_super, k=1, which='SA')
  78.     Energy = E[0]
  79.     EnergyPerBond = (Energy - LastEnergy) / 2;
  80.  
  81.     ## form the reduced density matrix by reshaping Psi into a matrix
  82.  
  83.     Dim = BlockH.shape[0]
  84.     PsiMatrix = np.mat(np.reshape(Psi, [Dim, Dim]))
  85.     Rho = PsiMatrix * PsiMatrix.H
  86.    
  87.     ## Diagonalize the density matrix
  88.     ## The eigenvalues are arranged in ascending order
  89.     D, V = np.linalg.eigh(Rho)
  90.  
  91.     ## Construct the truncation operator, which is the projector
  92.     ## onto the m largest eigenvalues of Rho
  93.  
  94.     T = np.mat(V[:, max(0,Dim-m):Dim])
  95.     TruncationError = 1 - np.sum(D[max(0,Dim-m):Dim])
  96.  
  97.     print("{:6} {:6} {:16.8f} {:12.8f} {:12.8f} {:12.8f}"
  98.           .format(i, 4+i*2, Energy, EnergyPerBond,
  99.                   ExactEnergy-EnergyPerBond, TruncationError))
  100.  
  101.     ## Truncate the block operators
  102.  
  103.     BlockH = T.H * BlockH * T
  104.     BlockSz = T.H * BlockSz * T
  105.     BlockSp = T.H * BlockSp * T
  106.     BlockSm = T.H * BlockSm * T
  107.     BlockI = T.H * BlockI * T
  108.  
  109. # finish
  110. print("Finished")
Get Code
Edit - History - Print - Recent Changes - Search
Page last modified on June 24, 2020, at 01:24 AM