Constructing an MPO from a sum of operator chains

This notebook demonstrates how to construct an MPO 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 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 MPO 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 lattice sites
nsites = 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 MPO 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 MPO

mpo = chemtensor.construct_mpo_from_opchains("double complex", nsites, chains, opmap, coeffmap, qsite)

mpo.bond_dims

[1, 2, 3, 3, 2, 1]

Such an MPO can now be used to perform, e.g., a DMRG calculation. Note that in this example, the MPO does not correspond to a Hermitian operator since one of the coefficients is complex.

mpo.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, chain, opmap, coeffmap) for chain in chains])

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

0.0