ITensor 学习笔记 (一)

这里主要是学习一下知名张量网络库 ITensor 的使用方法。我将从环境搭建开始记录学习的全过程。参考资料:

ITensor 官方网站:http://itensor.org

ITensor GitHub: https://github.com/ITensor/ITensor

一个张量网络的 Wiki: https://tensornetwork.org

环境搭建

这里使用 Windows10 Subsystem Ubuntu 作为运行环境,需要安装 g++、liblapack-dev、libblas-dev。然后再 checkout ITensor 的源代码:

git clone https://github.com/ITensor/ITensor itensor

随后应当把 makfile 的配置文件 options.mk.sample 复制为 options.mk。在这里编译器用的是 g++,因此 CCCOM 配置的值为:

CCCOM=g++ -m64 -std=c++17 -fconcepts -fPIC

而运行环境为:

PLATFORM=lapack BLAS_LAPACK_LIBFLAGS=-lpthread -L/usr/lib -lblas -llapack  

在 subsys 中输入 make 命令编译文件,此时会产生库文件 lib/libitensor.a,而我们应当在 tutorial 目录下开展自己的项目,其中还有一些教程案例。这里我们就从教程文件开始学习,并对其中的关键点进行说明。

Project: 01_one_site

在 ITensor 中最基本的类型为 Index 而不是 Tensor,并且 Tensor 也是由 Index 构成的,因此

auto s = Index(2,"s");

实际上定义了一个二维 Index,其 tag 为 $s$,一个 Index 可以有多个 tag 并使用逗号分割。以上代码是一个最简单的 Index 实例化的方法,更多的方法可以参见 官方文档。官方文章中还说明了一个重要的函数,即 prime() 函数,它的作用是复制一个新的 Index,并将其 PrimeLevel 递增 1,这其实可以看作是一种名称,其作用会在后面遇到。

现在我们再来看看如何实例化一个 ITensor 张量对象:

auto Sx = ITensor(s,prime(s));

ITensor 类型的构造函数需要传入的参数为任意多的 Index 类型对象,因此上面的代码即定义了一个 Index 为 $s$ 和 $\prime{s}$ 的张量。接下来我们需要设置这个张量的值,可以使用如下语句:

Sx.set(s=1,prime(s)=2,+0.5);
Sx.set(s=2,prime(s)=1,+0.5); 

即相当于 Index $s$ 取 $1$、并且 Index prime($s$) 取 $2$ 时,对应的值为 $0.5$。如果我们使用 PrintData(Sx) 会得到输出:

(2,1) 0.5000000
(1,2) 0.5000000

Project: 02_two_site

这个例子中,定义了三种张量,分别是 Sz、Sp 和 Sm,并用它们构建了一个 Heisenberg Hamiltonian。关于里面的一些问题还可以参考 官方文档

海森堡模型

我们先来回顾一下海森堡模型,参考资料点这里,海森堡提出自旋与自旋之间可能有交互作用,而海森堡模型的哈密顿量为:

$$H=\frac{1}{2} \sum\limits_{i,j}J_{ij} \vec{S}_i \cdot \vec{S}_j$$

这是最一般的哈密顿量,$i$ 和 $j$ 表示晶格中格点的序号,而 $S \in {\sigma_x, \sigma_y, \sigma_z}$ 为自旋算符(泡利算子),$J_{ij}$ 为 exchange constant,表示的是能量的大小,在这里只考虑最临近的自旋才具有相互作用,在这里我们只考虑两种情况,$i=1, j=2$,并且 $J_{12} = J$,因此:

$$H=J[\vec{S}_1 \cdot \vec{S}_2] = J[S_1^xS_2^x+S_1^yS_2^y+S_1^zS_2^z]$$

并定义升降算符 $S_j^{\pm} = S_j^x \pm iS_j^y$,因此哈密顿量可以改写为:

$$H=J[\frac{1}{2}(S_1^+S_2^- + S_1^-S_2^+) + S_1^zS_2^z]$$

这个就是我们需要构造的海森堡模型的哈密顿量的最终形式。

使用 ITensor 构造

接下来我们尝试使用 ITensor 来构造这个哈密顿量,过程很简单,这里只把关键代码列出来:

// 定义 Index
auto s1 = Index(2,"s1");
auto s2 = Index(2,"s2");
// 定义单格点算子
auto Sz1 = makeSz(s1);
auto Sz2 = makeSz(s2);
auto Sp1 = makeSp(s1);
auto Sp2 = makeSp(s2);
auto Sm1 = makeSm(s1);
auto Sm2 = makeSm(s2);
// 定义两格点的海森堡哈密顿量
auto H = Sz1*Sz2 + 0.5*Sp1*Sm2 + 0.5*Sm1*Sp2;

其实 * 表示矩阵乘法。最关键的在于三个单格点算子的定义方法,其中 p 和 m 分别代表 plus 和 minus,即上升算符和下降算符:

ITensor makeSp(Index const& s)
    {
    auto Sp = ITensor(s,prime(s));
    Sp.set(s=2,prime(s)=1, 1);
    return Sp;
    }
ITensor makeSm(Index const& s)
    {
    auto Sm = ITensor(s,prime(s));
    Sm.set(s=1,prime(s)=2,1);
    return Sm;
    }
ITensor makeSz(Index const& s)
    {
    auto Sz = ITensor(s,prime(s));
    Sz.set(s=1,prime(s)=1,+0.5);
    Sz.set(s=2,prime(s)=2,-0.5);
    return Sz;
    }

这里肯能不太好理解,我们使用一个图来把这个关系说清楚就好理解了:

计算初始能量期望值

其实也就是求哈密顿量 $H$ 在本征态 $|\psi\rangle$ 上的投影值。这里可以直接给出代码:

auto initEn = elt(dag(prime(psi)) * H * psi);
printfln("\nInitial energy = %.10f",initEn);

这里 ITensor dag(Itensor itensor) 函数是返回与 itensor 厄米共轭的张量。而 ITensor::Real elt(ITensor itensor) 则是返回$\langle \psi |H|\psi \rangle$的期望值,其实上面的这种写法本来就是返回一个标量,而 elt() 函数的作用主要在于一个类型转换,即把 ITensor 类型转换成 ITensor::Real,而它本身即为 double 类型。

求哈密顿量的基态

这里的参考资料主要为 QCQI 的第 4.7.2 节。代码如下:

Real beta = 0.;
auto expH = expHermitian(H,-beta);

// Here we apply exp(-beta*H), normalize
// and unprime
auto psibeta = expH*psi;
psibeta.noPrime();
psibeta /= norm(psibeta);

auto En = elt(dag(prime(psibeta)) * H * psibeta);
printfln("At beta=%.1f, energy = %.10f",beta,En);

这段代码首先构造了一个酉算子 $U=\exp(-\beta H)$,并作用到量子态 $|\psi\rangle$ 上,然后进行归一化。这里需要注意的是 noPrime() 函数,功能是将 prime level 设置为 $0$。现在,我们可以传入不同的 beta 并且得到变换后的能量值,如果 $\beta$ 足够大,对应的态即为基态。

另外官方文档中有 相关的例子

随后我们使用 SVD 进行分解,用法参考 官方文档

auto [U,D,V] = svd(psibeta,{s1});

用法很简单,需要注意的是,第二个参数 {s1} 代表把 s1 看作是行。另外,这里的 svd 函数还可以附带更多的参数,例如:

auto [U,S,V] = svd(psibeta, {s1}, {"Cutoff=",1E-2,"MaxDim=",1});

其中 Cutoff 表示去掉奇异值在其之下的项,而 MaxDim 表示最多保留的奇异值的数量。

ITensor 的 MPS

先来看看如何在 ITensor 中创建 MPS,可以参考 官方文档,以及 另一份官方文档。代码如下:

SpinHalf sites(5);
auto psi = MPS(sites);

这里使用 MPS 函数是生成随机的乘积态,我们可以使用 psi.A(…) 函数以只读形式取出其中的某个张量,我们取出第一个和第二个张量并打印来看看:

psi.A(1) =
ITensor ord=2:
(2|id=170|n=1,Site,S=1/2) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(1|id=700|l=1,Link) <Out>
  1: 1 QN()
{Zero / Not yet allocated}

psi.A(2) =
ITensor ord=3:
(1|id=700|l=1,Link) <In>
  1: 1 QN()
(2|id=281|n=2,Site,S=1/2) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(1|id=338|l=2,Link) <Out>
  1: 1 QN()
{Zero / Not yet allocated}

我们来进行简单的解释。 其中 (2|id=170|n=1,Site,S=1/2) 达标一个 Index,第一列代表维度,“n=1,Site,S=1/2” 是标签。其中 QN 表示量子数,参见 官方文档,其中 “Sz” 是名称,后面即为值,由于 $S=1/2$,因此 $QN=\pm 1$ 代表自旋为 $\pm 1/2$。显然,上面定义了一个链接维度为 1 的 MPS。

这里又一个问题指的注意,根据官方文档中的说法,使用构造函数 MPS() 就能生成随机的 MPS,但结果并非如此,例如上面的输出结果实际上只是生成了一个 MPS 的结构罢了,并且 “{Zero / Not yet allocated}” 也明确说明了这一点。而真正要生成一个具体的随机 MPS,还需要使用 randomMPS() 函数:

int N = 5;
SpinHalf sites(5, {"ConserveQNs=", false});
auto psi = randomMPS(sites);

我们看看输出结果:

ITensor ord=2: (2|id=565|n=1,Site,S=1/2) (1|id=817|l=1,Link)
{norm=0.99 (Dense Real)}
(1,1) 0.3385940
(2,1) 0.9334875

ITensor ord=3: (1|id=817|l=1,Link) (2|id=299|n=2,Site,S=1/2) (1|id=654|l=2,Link)
{norm=0.51 (Dense Real)}
(1,1,1) 0.0635747
(1,2,1) 0.5014503

...

这样就对了,只是,现在完全没有没有归一化啊。我们来看看 position(k) 函数的作用,它实际上就是移动正交中心 (orthogonality center) 到第 k 个张量上去,我们举个例子,例如我们运行 position(2),将会得到:

ITensor ord=2: (2|id=565|n=1,Site,S=1/2) (1|id=521|l=1,Link)
{norm=1.00 (Dense Real)}
(1,1) 0.3409816
(2,1) 0.9400700

ITensor ord=3: (2|id=299|n=2,Site,S=1/2) (1|id=442|l=2,Link) (1|id=521|l=1,Link)
{norm=0.05 (Dense Real)}
(1,1,1) 0.0062894
(2,1,1) 0.0496079

ITensor ord=3: (2|id=279|n=3,Site,S=1/2) (1|id=742|l=3,Link) (1|id=442|l=2,Link)
{norm=1.00 (Dense Real)}
(1,1,1) 0.9466765
(2,1,1) 0.3221856

...

关于正交中心和正则化的知识可以参考 这篇论文这篇论文 以及 这篇论文

Leave a Reply

Your email address will not be published. Required fields are marked *