Constructing a TTNO from a sum of operator chains

This notebook demonstrates how to construct a TTNO (tree tensor network operator) from a sum of operator chains of the form

\[\sum_j c_j \cdot \mathrm{op}_{j,0} \otimes \mathrm{op}_{j,1} \otimes \cdots \otimes \mathrm{op}_{j,n-1},\]

where \(n\) is the number of physical lattice sites, each \(\mathrm{op}_{j,\ell}\) is a local operator acting on site \(\ell\) (numbering starting from zero) and \(c_j\) is a coefficient. Internally, chemtensor optimizes the virtual bond dimensions of the TTNO based on bipartite graph theory. In this basic example, the local operators are the Pauli matrices and the identity.

import numpy as np
import chemtensor

# maximum number of OpenMP threads (0 indicates that OpenMP is not available)
chemtensor.get_max_openmp_threads()

16

# number of physical lattice sites
nsites_physical = 5

# physical quantum numbers at each site (all zero in this example)
qsite = [0, 0]

Define operator chains

# Pauli matrices
sigma_x = np.array([[0.,  1.], [1.,  0.]])
sigma_y = np.array([[0., -1j], [1j,  0.]])
sigma_z = np.array([[1.,  0.], [0., -1.]])

# operator map; the operator identifiers (OIDs) are the indices for this lookup-table, e.g., OID 2 refers to Pauli-Y
opmap = [np.identity(2), sigma_x, sigma_y, sigma_z]

# coefficient map; first two entries must always be 0 and 1;
# the coefficient identifiers (CIDs) are the indices for this lookup-table
coeffmap = [0, 1, 0.8, 0.4 - 1.5j, -0.7]

# operator chains, containing operator identifiers (OIDs) and coefficient identifiers (CIDs)
# referencing 'opmap' and 'coeffmap', respectively
chains = [
    #                   OIDs          qnumbers       CID start site
    chemtensor.OpChain([1, 0, 2],    [0, 0, 0, 0],    4, istart=2),  # coeffmap[4] * X_2 I_3 Y_4
    chemtensor.OpChain([1, 3, 0, 2], [0, 0, 0, 0, 0], 1, istart=0),  # coeffmap[1] * X_0 Z_1 I_2 Y_3
    chemtensor.OpChain([3, 3],       [0, 0, 0],       3, istart=1),  # coeffmap[3] * Z_1 Z_2
    chemtensor.OpChain([2],          [0, 0],          2, istart=4),  # coeffmap[2] * Y_4
]

In general, the integer quantum numbers are interleaved with the local operators to implement abelian symmetries (like particle number conservation). In practice, the quantum numbers endow the TTNO tensors with a sparsity pattern, handled internally by chemtensor. In this basic example, we do not use symmetries for simplicity and set all quantum numbers to zero.

Construct the TTNO

Tree topology:

#     0     4
#      ╲   ╱
#       ╲ ╱
#        5
#        │
#        │
#  1 ─── 6 ─── 3
#        │
#        │
#        2

tree_neighbors = [
    [5],           # neighbors of site 0
    [6],           # neighbors of site 1
    [6],           # neighbors of site 2
    [6],           # neighbors of site 3
    [5],           # neighbors of site 4
    [0, 4, 6],     # neighbors of site 5
    [1, 2, 3, 5],  # neighbors of site 6
]

Sites 0, …, 4 are physical sites, and the remaining sites 5, 6 are branching sites (with dummy physical legs of dimension 1).

ttno = chemtensor.construct_ttno_from_opchains("double complex", nsites_physical, tree_neighbors, chains, opmap, coeffmap, qsite)

# number of branching sites (in this example, sites 5 and 6)
ttno.nsites_branching

2

# show bond dimension between sites 1 and 6
ttno.bond_dim(1, 6)

2

# show bond dimension between sites 5 and 6
ttno.bond_dim(5, 6)

3

# get the list of TTNO tensors
ttno_tensors = ttno.a

# show one of the TTNO tensors
ttno_tensors[5]

array([[[[[1.+0.j, 0.+0.j, 0.+0.j]]],


        [[[0.+0.j, 1.+0.j, 0.+0.j]]]],



       [[[[0.+0.j, 0.+0.j, 0.+0.j]]],


        [[[0.+0.j, 0.+0.j, 1.+0.j]]]]])

# the two dimensions of 1 are the (dummy) physical dimensions, and the others the virtual bond dimensions (ordered by site index)
ttno_tensors[5].shape

(2, 2, 1, 1, 3)

# overall matrix representation of the TTNO
ttno.to_matrix()

array([[0.4-1.5j, 0. -0.8j, 0. +0.j , ..., 0. +0.j , 0. +0.j , 0. +0.j ],
       [0. +0.8j, 0.4-1.5j, 0. +0.j , ..., 0. +0.j , 0. +0.j , 0. +0.j ],
       [0. +0.j , 0. +0.j , 0.4-1.5j, ..., 0. +0.j , 0. +0.j , 0. +0.j ],
       ...,
       [0. +0.j , 0. +0.j , 0. +0.j , ..., 0.4-1.5j, 0. +0.j , 0. +0.j ],
       [0. +0.j , 0. +0.j , 0. +0.j , ..., 0. +0.j , 0.4-1.5j, 0. -0.8j],
       [0. +0.j , 0. +0.j , 0. +0.j , ..., 0. +0.j , 0. +0.8j, 0.4-1.5j]])

Reference calculation, as consistency check

def pad_op_chain(chain: chemtensor.OpChain, new_length: int):
    """
    Construct a new OpChain with identities padded on the left and right.
    """
    npad_right = new_length - chain.length - chain.istart
    assert npad_right >= 0
    # OID 0 always represents the local identity operation
    return chemtensor.OpChain(chain.istart*[0] + chain.oids  + npad_right*[0],
                              chain.istart*[0] + chain.qnums + npad_right*[0],
                              chain.cid, istart=0)

def op_chain_to_matrix(nsites: int, chain: chemtensor.OpChain, opmap, coeffmap):
    """
    Represent the logical operation of the operator chain as a dense matrix.
    """
    if chain.istart > 0 or chain.length < nsites:
        chain = pad_op_chain(chain, nsites)
    assert chain.istart == 0
    assert chain.length == nsites
    mat = coeffmap[chain.cid] * np.identity(1)
    for oid in chain.oids:
        mat = np.kron(mat, opmap[oid])
    return mat

# reference matrix representation
mat_ref = sum([op_chain_to_matrix(nsites_physical, chain, opmap, coeffmap) for chain in chains])

# compare (difference should be zero)
np.linalg.norm(ttno.to_matrix() - mat_ref)

0.0