Advanced Computing Platform for Theoretical Physics

Commit 7a91c2ad authored by Pan Zhang's avatar Pan Zhang
Browse files

CTM, iPEPS for 2D Ising

parent 71b5a85c
%% Cell type:markdown id: tags:
# iPEPS contraction for computing $\ln Z$ of the $2-$D ferromagnetic Ising model
%% Cell type:markdown id: tags:
## $Z$ as contraction of PEPS
The partition function of $2-$D Ising model can be written as
$$Z=\sum_\mathbf{s}\prod_{ij}e^{\beta J_{ij}s_is_j}=\sum_\mathbf{s}\prod_i^\circ\sqrt{A_i}.$$ Where $\prod_i^\circ$ denotes tensor contraction of all the $A_i$ tensors, each of which is a tensor given by $$A_i = I_i\prod_{j\in\partial i}\sqrt{B_{ij}},$$ where $I$ is a identity tensor with order equals to the degree of node $i$, and $B_{ij}$ defines a Boltzmann matrix with $$B_{ij}=\left[\begin{array}{l}e^{\beta J_{ij}}& e^{-\beta J_{ij}} \\ e^{-\beta J_{ij}}& e^{\beta J_{ij}}\end{array}\right].$$ $\sqrt{B}$ indicates square root of matrix $B$, yileding $\sqrt{B}\times \sqrt{B} = B$. For the ferromagnetic model, every $B$ matrix and $A$ tensor are identical, the system has a rotational invariance.
%% Cell type:code id: tags:
``` python
import torch
from scipy.linalg import sqrtm
import numpy as np
import sys,math
mydevice = torch.device('cpu')
mydtype=torch.float64
beta=0.5
def get_A(beta=0.5):
b=np.exp(beta)
B=torch.tensor(sqrtm(np.array([[b,1/b],[1/b,b]])),dtype=mydtype,device=mydevice) # the Boltzmann matrix
I=torch.zeros([2,2,2,2],dtype=mydtype) # the identity matrix defined on each site
I[0,0,0,0]=1
I[1,1,1,1]=1
Im=I.clone()
Im[1,1,1,1]=-1
A=torch.einsum("ij,ab,cd,xy,jbdy->iacx",B,B,B,B,I) #
Am= torch.einsum("ij,ab,cd,xy,jbdy->iacx",B,B,B,B,Im) # this will be used for computing magnetization
return A,Am
```
%% Cell type:markdown id: tags:
## Onsager solution
What we are going to compare with is the Onsager solution, with $\ln Z$, critical temperature $\beta_c$, and the spontaneous magnetization are given as
\begin{align} \ln Z&=\ln 2 +\frac{1}{8\pi^2}\int_0^{2\pi}d\theta_1\int_0^{2\pi}d\theta_2\ln\left[ \cosh^2(2\beta)-\sinh(2\beta)\cos(\theta_1)-\sinh(2\beta)\cos(\theta_2) \right]\\
\beta_c&=\frac{\ln (1+\sqrt{2})}{2},\\
m_{\beta \gt \beta_c}&=\left[ 1-\sinh^{-4}(2\beta) \right]^{\frac{1}{8}}.
\end{align}
%% Cell type:code id: tags:
``` python
def m_exact(beta):
if(beta>0.5*math.log(1+math.sqrt(2))):
return math.pow((1- math.pow(math.sinh(2*beta),-4)),1/8 )
else:
return 0
```
%% Cell type:markdown id: tags:
## CTMRG for rotational and translational symmetric iPEPS
The Corner Transfer Matrix RG method for contracting iPEPS takes care of a corner matrix $C$ and an edge matrix $E$, with $Z$ finally expressed as
$$ Z=\begin{array}{c}C&-&E&-&C\\|&&|&&|\\E&-&A&-&E\\|&&|&&|\\C&-&E&-&C\end{array}$$
The $C$ is updated by
$$ \begin{array}{c}C&-&E&-\\|&&|\\E&-&A&-\\|&&|&\\\end{array}\,\,\,\,\,\longrightarrow\begin{array}{c}\widehat C&=\\\|&\\ \end{array}\,\,\,\,\,\longrightarrow \begin{array}{c}C&-&U&=\\|&&&\\U^{\dagger}&&&\\\|&&& \end{array}$$
And we can see that after one step of update, $C$ is actually the diagonal matrix containing truncated eigenvalues of $\widehat C$, and $U$ contains conrreponding eigenvectors.
The edge tensor $E$ is updated as
$$ \begin{array}{c}|&&|\\E&-&A&-\\|&&| \end{array}\,\,\,\,\,\longrightarrow \begin{array}{c} \|\\ \widehat E&-\\ \| \end{array} \,\,\,\,\,\longrightarrow \,\,\,\,\, \begin{array}{c} \|\\ U \\ |\\U^{\dagger}\\ \| \\ \widehat E&-\\\|\\U\\|\\U^{\dagger}\\ \|\end{array} \longrightarrow \,\,\,\,\, \begin{array}{c} \|\\ U \\ |\\ E&-\\|\\U^{\dagger}\\\|\end{array}$$
For the detailed einsum rules for updates corner matrix $C$, I have used in the code that
$$ \begin{array}{c}C&-_j&E&-_k\\|_i&&|_l\\E&-_n&A&-_a\\|_m&&|_b&\\\end{array}$$
For update of $E$,
$$ \begin{array}{c}|_j&&|_l\\E&-_k&A&-_m\\|_i&&|_n \end{array}$$
And for finall contractions
$$ Z=\begin{array}{c}C&-_b&E&-_d&C\\|_a&&|_c&&|_e\\E&-_g&A&-_i&E\\|_f&&|_h&&|_j\\C&-_k&E&-_l&C\end{array}$$
Also be very careful about the order of the original indices of bonds when two bonds are merged into a singer boarder bond.
%% Cell type:code id: tags:
``` python
beta=0.48
A,Am = get_A(beta)
iter_max = 100
tol = 1.0e-10
Dmax = 10
C = B@torch.eye(2,dtype=mydtype)@B # initial C
I3=torch.zeros([2,2,2],dtype=mydtype)
I3[0,0,0]=1
I3[1,1,1]=1
E = torch.einsum("ij,ab,cd,jbd->iac",B,B,B,I3) # initial E
C = torch.randn(2,2,dtype=mydtype)
E = torch.randn(2,2,2,dtype=mydtype)
err=1.0
for iter in range(iter_max):
#update C
Chat = torch.einsum("ij,jkl,min,nlab->bmak",C,E,E,A).contiguous().view([E.shape[0]*A.shape[3],E.shape[1]*A.shape[2]])
[U,s,V]=torch.svd(Chat)
if(s.numel()>Dmax):
s=s[:Dmax]
U=U[:,:Dmax]
Cnew = torch.diag(s)
Cnew = Cnew/Cnew.norm()
if(C.shape == Cnew.shape):
err=(Cnew-C).norm()
#sys.stdout.write("#%d err=%.3g \n"%(iter,err))
if(err<tol):
#print("\nC matrix converged\n")
break
C=Cnew.clone()
E = torch.einsum("ijk,klmn->niljm",E,A).contiguous().view(E.shape[0]*A.shape[3],E.shape[1]*A.shape[1],A.shape[2])
E = torch.einsum("ijk,jl,im->mlk",E,U,U)
E = E/E.norm()
Z=torch.einsum("ab,bdc,de,fag,hgci,eji,kf,lkh,jl",C,E,C,E,A,E,C,E,C)
Zm=torch.einsum("ab,bdc,de,fag,hgci,eji,kf,lkh,jl",C,E,C,E,Am,E,C,E,C)
print("beta=%g,beta_C=%g"%(beta,0.5*math.log(1+math.sqrt(2))))
print("m_ctm = %g"%(Zm/Z).item(),"m_exact=%g"% m_exact(beta))
```
%%%% Output: stream
beta=0.48,beta_C=0.440687
m_ctm = -0.877523 m_exact=0.877523
%% Cell type:code id: tags:
``` python
```
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment