Let us try a simple idea, using the recursion scheme described above. At every step in the recursion, we add one spin on the right, and our basis dimension grows by a factor $2$. At some point during this recursion, the matrix will be too large to deal with. So let us fix a maximum number of states that we want to keep, $m$. At certain point during the process, the basis dimension will become larger than $m$. It is here that we start applying the truncation rule: diagonalize the Hamiltonian matrix exactly, and keep only the $m$ states with *lowest* eigenvalues(see Fig. [fig:nrg]).
As the system grows, the basis of the left block changes as we rotate to the new basis of eigenstates of the Hamiltonian. This is done by using a unitary transformation $U$. This matrix $U$ is nothing else but the matrix with the eigenstates ordered in columns. Therefore, adding a spin to the block now involves two steps: (i) we need to build the ‘tilde’ operators as before, and (ii) rotate the Hamiltonian matrix and the tilde operators to the new basis.
Let us assume that our old block before adding a site has a basis $\{|\alpha_{i-1}\rangle\}$, of dimension $D_{i-1}$, and the site has a basis $\{|s_{i}\rangle\}$ of dimension $d$. The new block basis $\{|\alpha_{i-1},s_{i}\rangle\}$ has dimension $d\times D_{i-1}$, such that we can easily diagonalize it to obtain all the eigenvalues and corresponding eigenvectors $\{|\alpha_{i+1}\rangle\}$. We build the matrix $U$ as the $D_{i-1} \times D_{i}$ unitary matrix with the $D_{i}=m$ eigenvectors with largest eigenvalues in the columns: $$U_{\alpha_{i-1}s_i,\alpha_i} = \langle \alpha_{i-1}s_i|\alpha_i\rangle.$$ Before the rotation, the operators had matrix elements: $$\tilde{O}_{\alpha_{i-1}s_i,\alpha'_{i-1}s'_i} = \langle \alpha_{i-1}s_i |\hat{O}|\alpha'_{i-1}s'_i\rangle.$$ We can now rotate all the tilde operators to the new basis as: $$\begin{aligned} \tilde{O}_{\alpha_i,\alpha'_i} & = & \langle \alpha_i |\hat{O}|\alpha'_i\rangle = \sum_{\alpha_{i-1},s_i} \sum_{\alpha'_{i-1}s'_i} \langle \alpha_i|\alpha_{i-1}s\rangle\langle \alpha_{i-1} s_i|\hat{O}|\alpha'_{i-1}s'_i\rangle\langle \alpha'_{i-1}s'|\alpha_i\rangle \nonumber \\ & = & \sum_{\alpha_{i-1},s_i} \sum_{\alpha'_{i-1}s'_i} (U^\dagger)_{\alpha_i,\alpha_{i-1}s_i} \tilde{O}_{\alpha_i \alpha'_i} U_{\alpha'_{i-1}s'_i,\alpha'_i}\end{aligned}$$ where the new matrices will have dimensions $m \times m$. we can now use these matrices to continue to block-growing process by adding another site. This can be repeated until the energy per site converges, or until we reach a desired system size.
It may seem that this new basis would be a natural choice is we assume that the physics of the problem is described by different manifolds with different energy scales. If we keep the lowest energy states and we get rid of the high energy states we can expect to get the low energy physics right. This in fact is the case in problems such as the Kondo and Anderson impurity problems@hewson. However, in strongly correlated, many-body problems such as the Heisenberg chain, this scheme performs poorly.
The problem was solved by Steve White by using what he called the ‘density matrix truncation’. He realized (without knowing at the time) that instead of getting rid of high energy states, one has to redistribute the ‘entanglement’ and minimize the loss of information. However, the way he formulated the problem did not incorporate the idea of entanglement, a concept that entered the picture much later after quantum information ideas were used to understand why and when the DMRG actually works. Before introducing these ideas, we shall describe the original formulation of the density matrix truncation@dmrg.
In order to introduce this new concept, we are going to use a new scheme: We are going to use two blocks instead of one, a left block, and a right block, as shown in Fig.[fig:two_blocks]. We are going to grow both blocks simultaneously using the same procedure outlined previously: at every step we add one site at the right of the left block, and one site to the left of the right block. The ground state can then be written as: $$|\Psi\rangle = \sum_{i,j} \Psi_{ij} |i\rangle |j\rangle, $$ where the sum runs over all the states of the left block $|i\rangle$ and right block $|j\rangle$, with the corresponding coefficients $\Psi_{ij}$.
Now the idea is as follows: once we reach the desired basis dimension $m$, we shall rotate the left block to a new basis $|i\rangle \rightarrow |\alpha\rangle$. We want to pick these states $|\alpha\rangle$ in such a way that when we truncate the basis, the “distance” between the original ground state $|\Psi\rangle$, and the new, truncated, variational approximation $|\tilde{\Psi}\rangle$, is minimized: $$S=\left||\Psi\rangle - |\tilde{\Psi}\rangle \right|^2,$$ where $$|\tilde{\Psi} \rangle = \sum_{\alpha=1}^m\sum_{j} \Psi_{\alpha j} |\alpha\rangle|j\rangle.$$
We are going to anticipate the solution: pick the basis $|\alpha\rangle$ given by the $m$ eigenvectors of the reduced density matrix of the left block with the $m$ largest eigenvalues. In order to justify this result, we first need to introduce some important concepts.
Imagine that we have a bipartite system, composed by subsystem $A$ and subsystem $B$, as shown in Fig.[fig:partition]. The Hilbert space of the system $A+B$ will be given by the tensor product of the Hilbert spaces of the two subsystems: $H_{A+B}=H_A \otimes H_B$, and will have dimension $D_{A+B} = D_A\times D_B$. Assume that the state of our system is described by a normalized wave-function $|\Psi\rangle$ that has support on $H_{A+B}$. We define the reduced density matrix of subsystem $A$ as $$\hat{\rho}_A = \mathrm{Tr}_B |\Psi\rangle\langle\Psi|.$$ Its corresponding matrix form is $${\rho_A}_{ii'} = \langle i|\hat{\rho}_A|i'\rangle = \sum_j \langle ij|\Psi\rangle\langle \Psi|i'j\rangle = \sum_{j} \Psi_{ij}\Psi^*_{i'j}$$ This operator describes the density matrix of a mixed state, in which the system $A$ is in contact with a bath or environment $B$. This is the price we have to pay for our ignorance of subsystem $B$.
The reduced density matrix has some nice properties:
It is Hermitian (or symmetric in case of real matrices). This means that its eigenvalues are real.
Its eigenvalues are non-negative.
The trace equals to unity: $\mathrm{Tr \rho_A} = 1$.
Its eigenvectors $|\alpha\rangle$ and eigenvalues $\omega_\alpha$ form an orthonormal basis.
This means that the reduced density matrix can be re-defined in the new eigenvector basis as: $$\hat{\rho}_A = \sum_\alpha \omega_\alpha |\alpha\rangle\langle\alpha|;$$ with $\omega_\alpha \ge 0$ and $\sum_\alpha {\omega_\alpha} = 1$.
These same considerations are valid for the block $B$.
Consider an arbitrary matrix $\Psi$ of dimensions $D_A \times D_B$. One can prove that $\Psi$ can de factorized as $$\Psi = UDV^{\dagger},$$ where $U$ is a ($D_A \times D_B$) unitary matrix, $V$ is a $(D_B \times D_B)$ unitary matrix, and $D$ is a $(D_B \times D_B$ diagonal matrix with real non-negative numbers along the diagonal, and zeroes elsewhere. Since $U$ and $V$ are unitary, they satisfy: $$\begin{aligned} UU^\dagger = {1}; \nonumber \\ VV^\dagger = {1}. \nonumber \end{aligned}$$ Their columns are orthonormal vectors, so $U$ and $V$ can be regarded as rotation matrices. The diagonal matrix elements $\lambda_\alpha$ of $D$ are known as the “singular values” of $\Psi$.
Let us apply the SVD to our quantum wave-function $|\Psi\rangle$ ([psi]), and for illustration, let us assume that $D_B \le D_A$. The coefficients $\Psi_{ij}$ define a matrix $\Psi$. After a SVD, they can be re-written as: $$\Psi_{ij} = \sum_\alpha^{D_B} U_{i\alpha}\lambda_\alpha(V^{\dagger})_{\alpha j} = \sum_\alpha^{D_B} U_{i\alpha}\lambda_\alpha V^{*}_{\alpha j}.$$ The wave-function can now be expressed as: $$\begin{aligned} |\Psi\rangle & = & \sum_{i}^{D_A} \sum_{j}^{D_B} \sum_\alpha^{D_B} U_{i\alpha}\lambda_\alpha V^*_{\alpha j}|i\rangle|j\rangle \nonumber \\ & = & \sum_\alpha^{D_B} \left(\sum_i^{D_A} U_{i\alpha}|i\rangle\right)\lambda_\alpha\left(\sum_j^{D_B} V^*_{\alpha j}|j\rangle\right) \nonumber \\ & = & \sum_\alpha^{D_B} \lambda_\alpha |\alpha\rangle_A|\alpha\rangle_B, \nonumber \end{aligned}$$ where we have defined the states $|\alpha\rangle_A = \sum_i U_{i\alpha}|i\rangle$ and $|\alpha\rangle_B = \sum_j V^*_{\alpha j}|j\rangle$. Due to the properties of $U$ and $V$, these states define a new orthogonal basis. This final expression is known as the “Schmidt decomposition” of the state $\Psi$, and the bases ${|\alpha\rangle}$ as the “Schmidt bases”..
In general, we have that the state $\Psi$ can be written in the new basis as: $$|\Psi\rangle = \sum_{\alpha}^r \lambda_\alpha |\alpha\rangle_A |\alpha\rangle_B; \,\,\, r=\mathrm{min}(D_A,D_B).$$
In the Schmidt basis, the reduced density matrices for the subsystems $A$ and $B$ are $$\rho_A = \mathrm{Tr}|\Psi\rangle\langle\Psi| = \sum_\alpha \lambda_\alpha^2|\alpha\rangle_{A}\ _A\langle \alpha|,$$ and $$\rho_B = \sum_\alpha \lambda_\alpha^2|\alpha\rangle_{B}\ _B\langle \alpha|$$
At this point, we realize some interesting observations:
The eigenvalues of the reduced density matrices are $\omega_\alpha = \lambda_\alpha^2$, the square of the singular values.
The two reduced density matrices share the spectrum.
The Schmidt bases are the eigenvectors of the reduced density matrices.
We here go back to the original problem of optimizing the wave-function in a reduced basis. In order to solve it, we are going to reformulate the question as: Given a matrix $\Psi$, what is the optimal matrix $\tilde{\Psi}$ with fixed rank $m$ that minimizes the Frobenius distance between the two matrices? It turns out, this is a well known problem called the “low ranking approximation”.
If we order the eigenvalues of the reduced density matrix in decreasing order $\omega_1,\omega_2,...\omega_m,...\omega_r$, it is straightforward to see that the Frobenius distance between the two matrices is given by $$S=\left|\Psi - \tilde{\Psi}\right|^2 = \sum_{m+1}^r \omega_i $$ This proves that the optimal basis is given by the eigenvectors of the reduced density matrix with the $m$ largest eigenvalues.
The above considerations allow us now to introduce the DMRG algorithm in a very natural way. We are going to present it in the traditional formulation, starting with the infinite-size algorithm, followed by the finite-size scheme.
The main idea behind the infinite-size algorithm consists in growing the left and right blocks by adding one site at a time. As we add sites, the basis of the blocks will grow, until we reach the desired maximum number of states $m$. At this point we need to start applying the density matrix truncation on both blocks. This process is repeated until we reach a desired system-size, or the error in the energy is below a pre-defined tolerance.
The algorithm illustrated in Fig.[fig:infinite] could be outlined as below:
Start growing the blocks by adding single-sites, as outlined in the exact diagonalization section. We assume that the Hilbert space for the single site has dimension $d$. When the size of the bocks become larger than $d \times m$, we start applying the density matrix truncation as follows:
In the early days of DMRG it was assumed that this scheme would lead to a good approximation of the system properties in the thermodynamic limit. Today we know that he best way to reach the thermodynamic limit is by using the finite-size algorithm on systems of fixed length, and doing a careful finite-size analysis of the results.
Let us now explain some of these steps in more detail.
Same as we did in the exact diagonalization section, we can add sites to the blocks by performing tensor products of the “tilde” operators on the block, and single-site operators.
Assume that we are in the $i^\mathrm{th}$ iteration of the algorithm, with our left and right blocks having length $i$. Let us label our $D_L$ basis states for the left block $\{|\alpha_{i}\rangle \}$, and our $d$ basis states for the single site that comes to the right $\{|s_{i+1}\rangle \}$ (See Fig.[fig:add_site]). When we add the site to the block, we obtain a new basis for the new combined block as $|\alpha_{i+1}\rangle = |\alpha_i\rangle \otimes |s_{i+1}\rangle$.
Let us assume for illustration purposes that we are dealing once more with the Heisenberg chain. All these ideas can be easily generalized to arbitrary models. Same as we did in the exact diagonalization section, we obtain the new Hamiltonian matrix for the combined block as: $$H_{L,i+1} = H_{L,i} \otimes {1}_2 + \tilde{S}^z_{L,i} \otimes S^z + \frac{1}{2} \left( \tilde{S}^+_{L,i} \otimes S^- + \tilde{S}^-_{L,i} \otimes S^+ \right).$$ In this expression, the “tilde” operators are in the $|\alpha_i\rangle$ basis, while the others are defined in the single-site basis.
A similar expression applies to the right block, which is obtained from the single site at position $i+2$, with basis $\{|s_{i+2}\rangle\}$ and dimension $d$, and the right block with basis $\{|\beta_{i+3}\rangle\}$ and dimension $D_R$: $$H_{R,i+2} = {1}_2 \otimes H_{R,i+3} + S^z \otimes \tilde{S}^z_{R,i+3} + \frac{1}{2} \left( {S^+} \otimes \tilde{S}^-_{R,i+3} + {S^-} \otimes \tilde{S}^+_{R,i+3} \right).$$
We now need to combine the left and right blocks to form the super-Hamiltonian: $$\hat{H}=\hat{H}_{L,i+1}+\hat{H}_{R,i+2} + S^z_{i+1} S^z_{i+2} + \frac{1}{2}\left( S^+_{i+1}S^-_{i+2} + S^-_{i+1}S^+_{i+2} \right)$$ where $\hat{H}_{L(R)}$ where obtained above, and only involve terms in the left (right) block. The single sites at positions $i+1$ and $i+2$ were absorbed by the left and right blocks, respectively, so in order to build the interactions, we have to rotate the corresponding operators to the new basis of the blocks. This is again done in the same spirit of the “tilde” transformation: $$\begin{aligned} H & = & H_{L,i+1} \otimes {1}_{D_R\times2} + {1}_{D_L \times 2} \otimes H_{R,i+2} \nonumber \\ & + & {1}_{D_L} \otimes S^z \otimes S^z \otimes {1}_{D_R} \nonumber \\ & + & \frac{1}{2} {1}_{D_L} \otimes S^+ \otimes S^- \otimes {1}_{D_R} \nonumber \\ & + & \frac{1}{2} {1}_{D_L} \otimes S^- \otimes S^+ \otimes {1}_{D_R} \nonumber\end{aligned}$$ or: $$\begin{aligned} H & = & H_{L,i+1} \otimes {1}_{D_R\times2} + {1}_{D_L \times 2} \otimes H_{R,i+2} \nonumber \\ & + & \tilde{S}^z_{L,i+1} \otimes \tilde{S}^z_{R,i+2} \nonumber \\ & + & \frac{1}{2} \tilde{S}^+_{L,i+1} \otimes \tilde{S}^-_{R,i+2} \nonumber \\ & + & \frac{1}{2} \tilde{S}^-_{L,i+1} \otimes \tilde{S}^+_{R,i+2} \nonumber \end{aligned}$$
The truncation process is similar to the one use in numerical renormalization group, but instead of using the matrix of eigenvectors of the Hamiltonian, we use the eigenvectors $\{|\alpha\rangle\},\{|\beta\rangle\}$ of the left and right reduced density matrix. Therefore, the new basis states for the left and right block are related to the states in the previous step as: $$\begin{aligned} |\alpha_{i+1}\rangle & = & \sum_{s_{i+1},\alpha_i}\langle \alpha_i s_{i+1}|\alpha_{i+1} \rangle |\alpha_i s_{i+1} \rangle = \sum_{s_{i+1},\alpha_i} (U_L)_{\alpha_i s_{i+1},\alpha_{i+1}} | \alpha_i s_{i+1} \rangle \nonumber \\ |\beta_{i+2}\rangle & = & \sum_{s_{i+2},\beta_{i+3}}\langle s_{i+2} \beta_{i+3}|\beta_{i+2} \rangle |s_{i+2} \beta_{i+3} \rangle = \sum_{s_{i+2},\beta_{i+3}} (U_R)_{s_{i+2} \beta_{i+3},\beta_{i+2}} | s_{i+2} \beta_{i+3} \rangle\end{aligned}$$ where $$(U_L)_{\alpha_i s_{i+1},\alpha_{i+1}} = \langle \alpha_i s_{i+1}|\alpha_{i+1} \rangle$$ and $$(U_R)_{s_{i+2} \beta_{i+3},\beta_{i+2}} = \langle s_{i+2} \beta_{i+3}|\beta_{i+2} \rangle.$$ If we keep only $m$ states, the matrices $U_{L(R)}$ will have dimensions $D_{L(R)} d\times m$. If the basis had already been truncated in the previous step, then $D_L(R)=m$.
We can now use these transformations to obtain the matrix elements for all the operators in the new truncated basis. For instance, an operator acting on a site inside the left block will be transformed as: $$\begin{aligned} (\tilde{O}_{L,i+1})_{\alpha_{i+1},\alpha'_{i+1}} & = & \langle \alpha_{i+1} |\hat{O}|\alpha'_{i+1}\rangle \nonumber \\ & = & \sum_{\alpha_i,s_{i+1}} \sum_{\alpha'_i,s'_{i+1}} \langle \alpha_{i+1}|\alpha_i s_{i+1}\rangle\langle \alpha_i s_{i+1}|\hat{O}|\alpha'_i s'_{i+1}\rangle\langle \alpha'_i s'_{i+1}|\alpha'_{i+1}\rangle \nonumber \\ & = & \sum_{\alpha_i s_{i+1}}\sum_{\alpha'_i s'_{i+1}} (U_L^\dagger)_{\alpha_{i+1},\alpha_is_{i+1}} (\tilde{O_{L,i}})_{\alpha_i s_{i+1},\alpha'_i s'_{i+1}} (U_L)_{\alpha'_is'_{i+1},\alpha'_{i+1}},\end{aligned}$$ and a similar expression can be obtained for operators in the right block.
As we mentioned before, the proper way to reach the thermodynamic limit with DMRG is by studying finite systems and performing a finite-size analysis. In order to study finite system, a generalization of the above ideas needs to be applied. The finite-size DMRG (illustrated in Fig.[fig:finite]) can be summarized as follows:
Once the desired system size is reached we start performing ``DMRG sweeps'', from right-to-left, and left-to-right to optimize the bases and improve accuracy. A left-to-right sweep is described as follows:
Perform a right-to-left sweep, by growing the right block one site at a time, and using the left block from the previous left-to-right sweep.
This sweeping process works in a similar fashion as a self-consistent loop, where we iteratively improve the solution. In fact, the DMRG can be formulated as a variational method, in which the variational parameters are continuously improved to minimize an energy functional. Intuitively a way to see it is by imagining a “demon” probing the environment around the block for the optimal states to improve the basis to represent the ground-state. These states are “absorbed” inside the block by the density matrix and its eigenvectors.
As described above, the shrinking block is replaced by the block from the previous sweep in the opposite direction. This means that all the information about the block and its operators needs to be stored, either in memory, or dumped on disk or stored in memory.
In [10]:
def psi_dot_psi(psi1, psi2):
x = 0.
for i in range(psi1.shape[0]):
for j in range(psi2.shape[1]):
x += psi1[i,j]*psi2[i,j]
return x
def lanczos(m, seed, maxiter, tol, use_seed = False, force_maxiter = False):
x1 = seed
x2 = seed
gs = seed
a = np.zeros(100)
b = np.zeros(100)
z = np.zeros((100,100))
lvectors = []
control_max = maxiter;
e0 = 9999
if(maxiter == -1):
force_maxiter = False
if(control_max == 0):
gs = 1
maxiter = 1
return(e0,gs)
x1[:,:] = 0
x2[:,:] = 0
gs[:,:] = 0
a[:] = 0.0
b[:] = 0.0
if(use_seed):
x1 = seed
else:
for i in range(x1.shape[0]):
for j in range(x1.shape[1]):
x1[i,j] = (2*np.random.random()-1.)
# x1[:,:] = 1
b[0] = psi_dot_psi(x1,x1)
b[0] = np.sqrt(b[0])
x1 = x1 / b[0]
x2[:] = 0
b[0] = 1.
e0 = 9999
nmax = min(99, maxiter)
for iter in range(1,nmax+1):
eini = e0
if(b[iter - 1] != 0.):
aux = x1
x1 = -b[iter-1] * x2
x2 = aux / b[iter-1]
aux = m.product(x2)
x1 = x1 + aux
a[iter] = psi_dot_psi(x1,x2)
x1 = x1 - x2*a[iter]
b[iter] = psi_dot_psi(x1,x1)
b[iter] = np.sqrt(b[iter])
lvectors.append(x2)
# print "Iter =",iter,a[iter],b[iter]
z.resize((iter,iter))
z[:,:] = 0
for i in range(0,iter-1):
z[i,i+1] = b[i+1]
z[i+1,i] = b[i+1]
z[i,i] = a[i+1]
z[iter-1,iter-1]=a[iter]
d, v = np.linalg.eig(z)
col = 0
n = 0
e0 = 9999
for e in d:
if(e < e0):
e0 = e
col = n
n+=1
e0 = d[col]
print "Iter = ",iter," Ener = ",e0
if((force_maxiter and iter >= control_max) or (iter >= gs.shape[0]*gs.shape[1] or iter == 99 or abs(b[iter]) < tol) or \
((not force_maxiter) and abs(eini-e0) <= tol)):
# converged
gs[:,:] = 0.
for n in range(0,iter):
gs += v[n,col]*lvectors[n]
print "E0 = ", e0
maxiter = iter
return(e0,gs) # We return with ground states energy
return(e0,gs)
In [15]:
import numpy as np
class Position:
LEFT, RIGHT = range(2)
class DMRGSystem(object):
def __init__(self, _nsites):
#Single site operators
self.nsites = _nsites
self.nstates = 2
self.dim_l = 0 # dimension of the left block
self.dim_r = 0 # dimension of the right block
self.left_size = 0 # number of sites in the left block
self.right_size = 0 # number of sites in the right block
self.sz0 = np.zeros(shape=(2,2)) # single site Sz
self.splus0 = np.zeros(shape=(2,2)) # single site S+
self.sz0[0,0] = -0.5
self.sz0[1,1] = 0.5
self.splus0[1,0] = 1.0
#Useful structures to store the matrices
self.HL = [] # left block Hamiltonian
self.HR = [] # right block Hamiltonian
self.szL = [] # left block Sz
self.szR = [] # right block Sz
self.splusL = [] # left block S+
self.splusR = [] # right block S+
zero_matrix = np.zeros(shape=(2,2))
for i in range(nsites):
self.HL.append(zero_matrix)
self.HR.append(zero_matrix)
self.szL.append(self.sz0)
self.szR.append(self.sz0)
self.splusL.append(self.splus0)
self.splusR.append(self.splus0)
self.psi = np.zeros(shape=(2,2)) # g.s. wave function
self.rho = np.zeros(shape=(2,2)) # density matrix
self.energy = 0.
self.error = 0.
#######################################
def BuildBlockLeft(self, iter):
self.left_size = iter
self.dim_l = self.HL[self.left_size-1].shape[0]
I_left = np.eye(self.dim_l)
I2 = np.eye(2)
# enlarge left block:
self.HL[self.left_size] = np.kron(self.HL[self.left_size-1],I2) + \
np.kron(self.szL[self.left_size-1],self.sz0) + \
0.5*np.kron(self.splusL[self.left_size-1],self.splus0.transpose()) + \
0.5*np.kron(self.splusL[self.left_size-1].transpose(),self.splus0)
self.splusL[self.left_size] = np.kron(I_left,self.splus0)
self.szL[self.left_size] = np.kron(I_left,self.sz0)
def BuildBlockRight(self, iter):
self.right_size = iter
self.dim_r = self.HR[self.right_size-1].shape[0]
I_right= np.eye(self.dim_r)
I2 = np.eye(2)
# enlarge right block:
self.HR[self.right_size] = np.kron(I2,self.HR[self.right_size-1]) + \
np.kron(self.sz0,self.szR[self.right_size-1]) + \
0.5* np.kron(self.splus0.transpose(),self.splusR[self.right_size-1]) + \
0.5* np.kron(self.splus0,self.splusR[self.right_size-1].transpose())
self.splusR[self.right_size] = np.kron(self.splus0,I_right)
self.szR[self.right_size] = np.kron(self.sz0,I_right)
def GroundState(self):
self.dim_l = self.HL[self.left_size].shape[0]
self.dim_r = self.HR[self.right_size].shape[0]
self.psi.resize((self.dim_l,self.dim_r))
maxiter = self.dim_l*self.dim_r
(self.energy, self.psi) = lanczos(self, self.psi, maxiter, 1.e-7)
def DensityMatrix(self, position):
# Calculate density matrix
if(position == Position.LEFT):
self.rho = np.dot(self.psi,self.psi.transpose())
else:
self.rho = np.dot(self.psi.transpose(),self.psi)
def Truncate(self, position, m):
# diagonalize rho
rho_eig, rho_evec = np.linalg.eigh(self.rho)
self.nstates = m
rho_evec = np.real(rho_evec)
rho_eig = np.real(rho_eig)
# calculate the truncation error for a given number of states m
# Reorder eigenvectors and trucate
index = np.argsort(rho_eig)
for e in index:
print "RHO EIGENVALUE ", rho_eig[e]
error = 0.
if (m < rho_eig.shape[0]):
for i in range(index.shape[0]-m):
error += rho_eig[index[i]]
print "Truncation error = ", error
aux = np.copy(rho_evec)
if (self.rho.shape[0] > m):
aux.resize((aux.shape[0],m))
n = 0
for i in range(index.shape[0]-1,index.shape[0]-1-m,-1):
aux[:,n]=rho_evec[:,index[i]]
n += 1
rho_evec = aux
# rho_evec = np.eye(self.rho.shape[0])
# perform transformation:
U = rho_evec.transpose()
if(position == Position.LEFT):
aux2 = np.dot(self.HL[self.left_size],rho_evec)
self.HL[self.left_size] = np.dot(U,aux2)
aux2 = np.dot(self.splusL[self.left_size],rho_evec)
self.splusL[self.left_size] = np.dot(U,aux2)
aux2 = np.dot(self.szL[self.left_size],rho_evec)
self.szL[self.left_size] = np.dot(U,aux2)
else:
aux2 = np.dot(self.HR[self.right_size],rho_evec)
self.HR[self.right_size] = np.dot(U,aux2)
aux2 = np.dot(self.splusR[self.right_size],rho_evec)
self.splusR[self.right_size] = np.dot(U,aux2)
aux2 = np.dot(self.szR[self.right_size],rho_evec)
self.szR[self.right_size] = np.dot(U,aux2)
def product(self, psi):
npsi = np.dot(self.HL[self.left_size],psi)
npsi += np.dot(psi,self.HR[self.right_size].transpose())
# Sz.Sz
tmat = np.dot(psi,self.szR[self.right_size].transpose())
npsi += np.dot(self.szL[self.left_size],tmat)
# S+.S-
tmat = np.dot(psi,self.splusR[self.right_size])*0.5
npsi += np.dot(self.splusL[self.left_size],tmat)
# S-.S+
tmat = np.dot(psi,self.splusR[self.right_size].transpose())*0.5
npsi += np.dot(self.splusL[self.left_size].transpose(),tmat)
return npsi
In [16]:
nsites = 10
n_states_to_keep = 10
n_sweeps = 4
S = DMRGSystem(nsites)
###############################################################################
for iter in range(1,nsites/2): # do infinite size dmrg for warmup
print "WARMUP ITERATION ", iter, S.dim_l, S.dim_r
# Create HL and HR by adding the single sites to the two blocks
S.BuildBlockLeft(iter)
S.BuildBlockRight(iter)
# find smallest eigenvalue and eigenvector
S.GroundState()
# Calculate density matrix
S.DensityMatrix(Position.LEFT)
# Truncate
S.Truncate(Position.LEFT,n_states_to_keep)
# Reflect
S.DensityMatrix(Position.RIGHT)
S.Truncate(Position.RIGHT,n_states_to_keep)
first_iter = nsites/2
for sweep in range(1,n_sweeps):
for iter in range(first_iter, nsites-3):
print "LEFT-TO-RIGHT ITERATION ", iter, S.dim_l, S.dim_r
# Create HL and HR by adding the single sites to the two blocks
S.BuildBlockLeft(iter)
S.BuildBlockRight(nsites-iter-2)
# find smallest eigenvalue and eigenvector
S.GroundState()
# Calculate density matrix
S.DensityMatrix(Position.LEFT)
# Truncate
S.Truncate(Position.LEFT,n_states_to_keep)
first_iter = 1;
for iter in range(first_iter, nsites-3):
print "RIGHT-TO-LEFT ITERATION ", iter, S.dim_l, S.dim_r
# Create HL and HR by adding the single sites to the two blocks
S.BuildBlockRight(iter);
S.BuildBlockLeft(nsites-iter-2)
# find smallest eigenvalue and eigenvector
S.GroundState();
# Calculate density matrix
S.DensityMatrix(Position.RIGHT)
# Truncate
S.Truncate(Position.RIGHT,n_states_to_keep)
In [ ]:
In [ ]:
In [ ]: