Usage examples

Holstein model

This example shows how to construct Hamiltonian of the Holstein model on a square lattice and to check that the total number of electrons is a conserved quantity of the model.

The Hamiltonian considered here is a sum of three terms.

• An electronic tight-binding model on a square lattice with only nearest-neighbour hopping allowed,

  \[\hat H_\text{e} = -t \sum_\sigma \sum_{\langle i,j\rangle} c^\dagger_{i,\sigma} c_{j,\sigma}.\]

• A harmonic oscillator at each lattice site (a localized phonon),

  \[\hat H_\text{ph} = \omega_0 \sum_i a^\dagger_i a_i.\]

• A coupling between the electrons and the phonons.

  \[\hat H_\text{e-ph} = g \sum_\sigma \sum_i n_{i,\sigma}(a^\dagger_i + a_i).\]

Instead of writing the sums over lattice sites explicitly, we call library functions tight_binding(), dispersion() and holstein_int(). We also make use of the NetworkX package to easily generate the adjacency matrix of the periodic square lattice.

#
# Holstein model on a square lattice with periodic boundary conditions.
#

from itertools import product
import numpy as np

from pycommute.expression import ExpressionR, conj, FERMION, BOSON
from pycommute.expression import n
from pycommute.models import tight_binding, dispersion, holstein_int

from networkx.generators.lattice import grid_2d_graph
from networkx.linalg.graphmatrix import adjacency_matrix

#
# Let us define Hamiltonian of an electronic tight-binding model
# on a square lattice.
#

# Number of lattice sites in each direction
# (the total number of sites is N*N).
N = 10

# Electron hopping constant - energy parameter of the TB model.
t = 2.0

# Use NetworkX to construct the periodic square lattice (graph)
lat = grid_2d_graph(N, N, periodic=True)

# Create lists of indices for electronic spin-up and spin-down operators.
# lat.nodes() returns a list of N^2 pairs of site indices (x, y).
indices_up = [(x, y, "up") for x, y in lat.nodes()]
indices_dn = [(x, y, "down") for x, y in lat.nodes()]

# A sum of tight-binding Hamiltonians for both spins.
# The hopping matrix passed to tight_binding() is proportional to the
# adjacency matrix of the lattice.
hopping_matrix = -t * adjacency_matrix(lat).todense()
H_e = tight_binding(hopping_matrix, indices=indices_up) \
    + tight_binding(hopping_matrix, indices=indices_dn)

#
# Hamiltonian of phonons localized at lattice sites.
#

# Frequency of the localized phonon.
w0 = 0.5

# A lists of indices for bosonic operators, simply the (x, y) pairs
indices_phonon = lat.nodes()

# All N^2 phonons have the same frequency
phonon_freqs = w0 * np.ones(N ** 2)
H_ph = dispersion(phonon_freqs, indices=indices_phonon, statistics=BOSON)

#
# Hamiltonian of electron-phonon coupling.
#

# Electron-phonon coupling constant.
g = 0.1

H_e_ph = holstein_int(g * np.ones(N ** 2),
                      indices_up=indices_up,
                      indices_dn=indices_dn,
                      indices_boson=indices_phonon)

# Complete Holstein Hamiltonian.
H_H = H_e + H_ph + H_e_ph

# Print H_H. There will be quite a lot of terms for the 100-site lattice!
print("H_H =", H_H)

# Check hermiticity of H_H.
print(r"H_H - H_H^\dagger =", H_H - conj(H_H))

# Check that H_H commutes with the total number of electrons N_e.
N_e = ExpressionR()
for spin in ("up", "down"):
    for x, y in product(range(N), range(N)):
        N_e += n(x, y, spin)

print("[H_H, N_e] =", H_H * N_e - N_e * H_H)

#
# Iteration interface
#

# Iterate over monomial-coefficient pairs in polynomial expression H_H
for monomial, coeff in H_H:
    print("Coefficient:", coeff)
    # Iterate over algebra generators (creation/annihilation operators)
    # in the monomial
    for generator in monomial:
        # Detect what algebra this generator belongs to
        if generator.algebra_id == FERMION:
            print("\tFermionic", end=' ')
        elif generator.algebra_id == BOSON:
            print("\tBosonic", end=' ')
        # Creation or annihilation operator?
        if generator.dagger:
            print("creation operator", end=' ')
        else:
            print("annihilation operator", end=' ')
        # Extract indices carried by the generator
        print("with indices", list(generator.indices))
        # N.B. generator.indices is an instance of a special tuple-like type
        # `pycommute.expression.Indices`. Types of its elements are restricted
        # to integer and string, and its comparison operators behave differently
        # from those of the Python tuple. Otherwise, it supports `len()`,
        # indexed element access and iteration protocol.

Spin-1/2 Heisenberg chain

The spin-1/2 Heisenberg chain is a textbook example of an integrable quantum system. Its Hamiltonian

\[\hat H = g \sum_i \mathbf{S}_i \cdot \mathbf{S}_{i+1}\]

conserves three projections of the total spin

\[\mathbf{S} = \sum_i \mathbf{S}_i\]

as well as a series of higher order charges \(Q_n\). Existence of these charges can be derived from the transfer matrix theory. Explicit expressions for \(Q_n\) were obtained in [GM94]. The following script constructs Hamiltonian of the Heisenberg chain with periodic boundary conditions (see heisenberg()) and checks that \([\hat H, \mathbf{S}] = 0\), \([\hat H, Q_n] = 0\) and \([Q_n, Q_m] = 0\) for \(m,n = 3,4,5\).

#
# Periodic spin-1/2 Heisenberg chain and its integrals of motion
#
# Expressions for the integrals of motion are taken from
#
#   "Quantum Integrals of Motion for the Heisenberg Spin Chain",
#   M. P. Grabowski and P. Mathieu
#   Mod. Phys. Lett. A, Vol. 09, No. 24, pp. 2197-2206 (1994),
#   https://doi.org/10.1142/S0217732394002057
#

from numpy import array, zeros, dot, cross
from pycommute.expression import S_x, S_y, S_z
from pycommute.models import heisenberg

# Number of spins in the chain
N = 20
# Heisenberg exchange constant
g = 2

# List of 3-component spin vectors {S_0, S_1, ..., S_{N-1}}
S = [array([S_x(i), S_y(i), S_z(i)]) for i in range(N)]

# Matrix of exchange constants between spins i and j
exchange_matrix = zeros((N, N))
# Set elements corresponding to the nearest neighbors to -g
# (index shift modulo N ensures periodic boundary conditions).
for i in range(N):
    exchange_matrix[i, (i + 1) % N] = -g

# Hamiltonian of the spin-1/2 Heisenberg chain.
H = heisenberg(exchange_matrix)

# Total spin of the chain.
S_tot = array(sum(S))

# All three components of S commute with the Hamiltonian.
print("[H, S_x] =", (H * S_tot[0] - S_tot[0] * H))
print("[H, S_y] =", (H * S_tot[1] - S_tot[1] * H))
print("[H, S_z] =", (H * S_tot[2] - S_tot[2] * H))

# Higher charge Q_3 (1st line of Eq. (10)).
Q3 = sum(dot(cross(S[i], S[(i + 1) % N]), S[(i + 2) % N]) for i in range(N))
print("[H, Q3] =", (H * Q3 - Q3 * H))

# Higher charge Q_4 (2nd line of Eq. (10)).
Q4 = sum(4.0 * dot(
    cross(
        cross(S[i], S[(i + 1) % N]), S[(i + 2) % N]
    ), S[(i + 3) % N]) for i in range(N))
Q4 += sum(dot(S[i], S[(i + 2) % N]) for i in range(N))
print("[H, Q4] =", (H * Q4 - Q4 * H))

# Higher charge Q_5 (3rd line of Eq. (10)).
Q5 = sum(4.0 * dot(
    cross(
        cross(
            cross(S[i], S[(i + 1) % N]),
            S[(i + 2) % N]
        ), S[(i + 3) % N]),
    S[(i + 4) % N]) for i in range(N))
Q5 += sum(dot(cross(S[i], S[(i + 2) % N]), S[(i + 3) % N]) for i in range(N))
Q5 += sum(dot(cross(S[i], S[(i + 1) % N]), S[(i + 3) % N]) for i in range(N))
print("[H, Q5] =", (H * Q5 - Q5 * H))

# Check that the higher charges pairwise commute.
print("[Q3, Q4] =", (Q3 * Q4 - Q4 * Q3))
print("[Q3, Q5] =", (Q3 * Q5 - Q5 * Q3))
print("[Q4, Q5] =", (Q4 * Q5 - Q5 * Q4))

Spectrum of Tavis-Cummings model

In this example we demonstrate how to use tools from the loperator module along with NumPy’s eigensolvers to diagonalize Hamiltonians of simple quantum models.

We will consider a two-atom Jaynes-Cummings (Tavis-Cummings) Hamiltonian [TC68], which models a system of two qubits coherently coupled via a bosonic degree of freedom (quantum oscillator),

\[\hat H = \epsilon (\hat S_{z,0} + \hat S_{z,1}) + \omega \hat a^\dagger \hat a + g_0 (\hat a^\dagger \hat S_{-,0} + \hat a \hat S_{+,0}) + g_1 (\hat a^\dagger \hat S_{-,1} + \hat a \hat S_{+,1}).\]

An expression of this form can be conveniently generated by a single call to jaynes_cummings().

Since the occupation number of the bosonic mode can formally grow to \(+\infty\), we have to artificially restrict the dimension of the bosonic Hilbert space in order to be able to construct a finite matrix representation of our problem.

Note

Performance of repeated calls to LOperatorR or to LOperatorC in order to construct a matrix representation of a linear operator is limited by Python method call overhead. For large-scale problems it is advised to call the utility function make_matrix().

#
# Diagonalization of a two-qubit Tavis-Cummings model.
#

import numpy as np

# Utility function used to construct the generalized Jaynes-Cummings model.
from pycommute.models import jaynes_cummings

# Hilbert spaces.
from pycommute.loperator import HilbertSpace, make_space_boson, make_space_spin

# Real-valued linear operator object.
from pycommute.loperator import LOperatorR

# Build the matrix form of a linear operator.
from pycommute.loperator import make_matrix

# Transition frequencies of qubits.
eps = np.array([1.0, 1.0])

# Oscillator frequency.
# (In a more general case of M oscillators, omega would be a length-M array).
omega = np.array([0.8])

# Qubit-oscillator coupling constants as a 2x1 array.
g = np.array([[0.5], [0.6]])

# Create the Tavis-Cummings Hamiltonian.
H = jaynes_cummings(eps, omega, g)

# Construct state space of our problem as a direct product of two
# two-dimensional Hilbert spaces (qubits) and one truncated bosonic Hilbert
# space.
# make_space_boson(4) returns the truncated bosonic space with allowed
# occupation numbers N = 0, 1, ..., (2^4-1).
hs = HilbertSpace([
    make_space_spin(1 / 2, 0),  # Qubit 1: spin-1/2, index 0
    make_space_spin(1 / 2, 1),  # Qubit 2: spin-1/2, index 1
    make_space_boson(4, 0)      # Oscillator, index 0
])

# Construct a linear operator corresponding to 'H' and acting in the Hilbert
# space 'hs'.
H_op = LOperatorR(H, hs)

#
# Prepare a matrix representation of 'H_op'
#

# Method I (manual).
H_mat1 = np.zeros((hs.dim, hs.dim))
for i in range(hs.dim):
    # A column vector psi = {0, 0, ..., 1, ..., 0}
    psi = np.zeros(hs.dim)
    psi[i] = 1.0
    # Act with H_op on psi and store the result the i-th column of H_mat
    H_mat1[:, i] = H_op * psi

# Method II (automatic and faster).
H_mat2 = make_matrix(H_op, hs)
print("Max difference between H_mat and H_mat2:",
      np.max(np.abs(H_mat1 - H_mat2)))

# Use NumPy to compute eigenvalues of H_mat
E = np.linalg.eigvals(H_mat1)

print("Energies:", np.sort(E))

Sectors of a multi-orbital interaction Hamiltonian

When implementing an exact diagonalization algorithm, it is usually advantageous to take an extra preparatory step and split the Hilbert space of the problem into subspaces (sectors) characterized by a set of quantum numbers. In many situations, it is possible to find the sectors even without knowing the quantities conserved by the Hamiltonian (see, for instance, [SKFP16]).

In the script shown below we construct an electron-electron interaction Hamiltonian of an atomic \(d\)-shell using Slater parametrization. Then, we employ the SpacePartition and BasisMapper utility classes (and optionally the make_matrix() function) to independently diagonalize the Hamiltonian within each sector.

#
# Sector-wise diagonalization of the Slater interaction Hamiltonian.
#

import numpy as np

# Utility function used to construct the Slater interaction Hamiltonian.
from pycommute.models import slater_int

# Hilbert space.
from pycommute.loperator import HilbertSpace

# Real-valued linear operator object.
from pycommute.loperator import LOperatorR
# Automatic space partition and basis state mapper
from pycommute.loperator import make_space_partition, BasisMapper
# Build the matrix form of a linear operator
from pycommute.loperator import make_matrix

# Values of Slater radial integrals F^0, F^2, F^4
F = np.array([4.0, 1.0, 0.2])

# Build an expression for the interaction Hamiltonian.
H = slater_int(F)

# Analyze structure of 'H' and construct a suitable Hilbert space.
hs = HilbertSpace(H)

# Construct a linear operator corresponding to H and acting in the Hilbert
# space 'hs'.
H_op = LOperatorR(H, hs)

# The first value returned by 'make_space_partition()' is a
# SpacePartition object. It represents a partition of the full Hilbert space
# into sectors, i.e. subspaces invariant under action of 'H_op'.
#
# The second returned value is a dictionary {(i, j): value} of all non-vanishing
# matrix elements of H_{ij}. By definition, all matrix elements of 'H_op'
# between different sectors vanish.
sp, matrix_elements = make_space_partition(H_op, hs)

# Print out information about the revealed sectors.
print("Dimension of full Hilbert space:", sp.dim)
print("Total number of sectors:", sp.n_subspaces)
for i in range(sp.dim):
    print("Many-body basis state %d belongs to sector %d" % (i, sp[i]))

# Compile lists of basis states spanning each sector.
sectors = sp.subspace_bases()

# Diagonalize 'H_op' within each sector
for n, sector in enumerate(sectors):
    sector_dim = len(sector)
    print("Diagonalizing sector %d of H_op" % n)
    print("Sector size is %d" % sector_dim)

    # A BasisMapper object translates indices of basis states from the full
    # Hilbert space to a given sector.
    basis_mapper = BasisMapper(sector)

    #
    # Prepare a matrix representation of 'H_op' within current sector.
    #

    # Method I (manual).
    H_mat1 = np.zeros((sector_dim, sector_dim))
    for i in range(sector_dim):
        # A column vector psi = {0, 0, ..., 1, ..., 0}
        psi = np.zeros(sector_dim)
        psi[i] = 1.0
        # This vector will receive the result of H_op * psi
        phi = np.zeros(sector_dim)

        # Since both 'psi' and 'phi' are 1D NumPy arrays with size of the chosen
        # sector, 'H_op' cannot act on them directly as it expects vectors of
        # size hs.dim.
        # Instead, we are going to use our basis mapper object to construct
        # special views of psi and phi.
        psi_view = basis_mapper(psi)
        phi_view = basis_mapper(phi)

        # Now, 'H_op' can act on the mapped views.
        H_op(psi_view, phi_view)

        # Store H_op * psi in the i-th column of 'H_mat1'.
        H_mat1[:, i] = phi

    # Method II (automatic and faster).
    H_mat2 = make_matrix(H_op, sector)
    print("Max difference between H_mat and H_mat2:",
          np.max(np.abs(H_mat1 - H_mat2)))

    # Use NumPy to compute eigenvalues of 'H_mat1'.
    E = np.linalg.eigvals(H_mat1)
    print("Energies:", np.sort(E))

N-fermion sector views

This example shows how to partially diagonalize the single band Fermi-Hubbard model on a large 2D square lattice with periodic boundary conditions.

\[\hat H = -t\sum_{\langle i,j \rangle, \sigma} c^\dagger_{i,\sigma} c_{j,\sigma} -\mu \sum_{i, \sigma} n_{i,\sigma} + U \sum_i n_{i,\uparrow} n_{i,\downarrow}.\]

The number of atoms in the considered lattice (cluster) is set to 16, which makes for an intractably large Hilbert space of dimension \(2^{32}\). Since the Hamiltonian of the Hubbard model preserves total numbers of spin-up and spin-down electrons, one can diagonalize the Hamiltonian within a single \(N\)-fermion sector or within a \((N_\uparrow, N_\downarrow)\)-multisector.

An \(N\)-fermion sector is a subspace of a full Hilbert space, which is spanned by all basis states with a fixed total occupation of fermionic degrees of freedom (FDOF). Similarly, a multisector is spanned by all basis states with a fixed occupation \(N_1\) of a subset of the FDOF \(\{S_1\}\), occupation \(N_2\) of another subset \(\{S_2\}\) and so on. There can be any number of pairs \((\{S_i\}, N_i)\) (sectors contributing to the multisector) as long as all the subsets \(\{S_i\}\) are disjoint.

In our example we consider a moderately sized sector with \(N = 2\) (\(\dim = {32 \choose 2} = 496\)) and a multisector with \(N_\uparrow = 1, N_\downarrow = 1\) (\(\dim = {16 \choose 1}{16 \choose 1} = 256\)).

Note

In general, the N-fermion (multi)sector views and functions do not require a purely fermionic system. The definitions of a sector and a multisector stand in presence of bosons or spins.

import numpy as np

# Utility functions used to construct the Fermi-Hubbard Hamiltonian.
from pycommute.models import tight_binding, dispersion, hubbard_int

# Hilbert space.
from pycommute.loperator import HilbertSpace

# Real-valued linear operator object.
from pycommute.loperator import LOperatorR
# Real-valued N-fermion (multi)sector views and related functions
from pycommute.loperator import (
    NFermionSectorViewR,
    NFermionMultiSectorViewR,
    n_fermion_sector_size,
    n_fermion_multisector_size
)

from networkx.generators.lattice import grid_2d_graph
from networkx.linalg.graphmatrix import adjacency_matrix

#
# Let us define Hamiltonian of a tight-binding model on a square lattice.
#

# Number of lattice sites in each direction
# (the total number of sites is Nx * Ny).
Nx = 4
Ny = 4

# Hopping constant
t = 0.5
# Chemical potential
mu = 1.0
# Coulomb repulsion
U = 2.0

# Use NetworkX to construct the periodic square lattice (graph)
lat = grid_2d_graph(Nx, Ny, periodic=True)

# Create lists of indices for spin-up and spin-down operators.
# lat.nodes() returns a list of Nx * Ny pairs of site indices (x, y).
indices_up = [(x, y, "up") for x, y in lat.nodes()]
indices_dn = [(x, y, "down") for x, y in lat.nodes()]

# A sum of tight-binding Hamiltonians for both spins.
# The hopping matrix passed to tight_binding() is proportional to the
# adjacency matrix of the lattice.
hopping_matrix = -t * adjacency_matrix(lat).todense()
H = tight_binding(hopping_matrix, indices=indices_up) \
    + tight_binding(hopping_matrix, indices=indices_dn)

# Add the chemical potential terms
H += dispersion(-mu * np.ones(len(indices_up)), indices=indices_up)
H += dispersion(-mu * np.ones(len(indices_dn)), indices=indices_dn)

# Add the Hubbard interaction term
H += hubbard_int(U * np.ones(len(indices_up)),
                 indices_up=indices_up,
                 indices_dn=indices_dn)

# Analyze structure of H and construct a suitable Hilbert space.
hs = HilbertSpace(H)
print("Full Hilbert space dimension:", hs.dim)

# Construct a linear operator corresponding to H and acting in the Hilbert
# space 'hs'.
H_op = LOperatorR(H, hs)

#
# Diagonalize the N = 2 sector of the model using 'NFermionSectorViewR'
#

N = 2

sector_size = n_fermion_sector_size(hs, N)
print("Size of the N = 2 sector:", sector_size)

# Prepare a matrix representation of H_op within the N = 2 sector.
H_mat = np.zeros((sector_size, sector_size))
for i in range(sector_size):
    # A column vector psi = {0, 0, ..., 1, ..., 0}
    psi = np.zeros(sector_size)
    psi[i] = 1.0
    # This vector will receive the result of H_op * psi
    phi = np.zeros(sector_size)

    # Since both psi and phi are 1D NumPy arrays with size of the chosen
    # sector, H_op cannot act on them directly as it expects vectors of
    # size hs.dim.
    # Instead, we are going to use real-valued 2-fermion sector views of
    # psi and phi.
    psi_view = NFermionSectorViewR(psi, hs, N)
    phi_view = NFermionSectorViewR(phi, hs, N)

    # Now, H_op can act on the mapped views.
    H_op(psi_view, phi_view)

    # Store H_op * psi in the i-th column of H_mat.
    H_mat[:, i] = phi

# Use NumPy to compute eigenvalues of H_mat.
E = np.linalg.eigvals(H_mat).real
print("10 lowest eigenvalues of the N = 2 sector:", np.sort(E)[:10])

#
# Diagonalize the N_up = 1, N_down = 1 multisector of the model using
# 'NFermionMultiSectorViewR'
#

N_up = 1
N_dn = 1

# Define sectors contributing to the multisector
sector_up = (indices_up, N_up)
sector_dn = (indices_dn, N_dn)

sectors = [sector_up, sector_dn]

multisector_size = n_fermion_multisector_size(hs, sectors)
print("Size of the N_up = 1, N_down = 1 multisector:", multisector_size)

# Prepare a matrix representation of H_op within the multisector.
H_mat = np.zeros((multisector_size, multisector_size))
for i in range(multisector_size):
    # A column vector psi = {0, 0, ..., 1, ..., 0}
    psi = np.zeros(multisector_size)
    psi[i] = 1.0
    # This vector will receive the result of H_op * psi
    phi = np.zeros(multisector_size)

    # Since both psi and phi are 1D NumPy arrays with size of the chosen
    # sector, H_op cannot act on them directly as it expects vectors of
    # size hs.dim.
    # Instead, we are going to use real-valued multisector views of psi and phi.
    psi_view = NFermionMultiSectorViewR(psi, hs, sectors)
    phi_view = NFermionMultiSectorViewR(phi, hs, sectors)

    # Now, H_op can act on the mapped views.
    H_op(psi_view, phi_view)

    # Store H_op * psi in the i-th column of H_mat.
    H_mat[:, i] = phi

# Use NumPy to compute eigenvalues of H_mat.
E = np.linalg.eigvals(H_mat).real
print("10 lowest eigenvalues of the multisector:", np.sort(E)[:10])

References

[GM94]

“Quantum Integrals of Motion for the Heisenberg Spin Chain”, M. P. Grabowski and P. Mathieu, Mod. Phys. Lett. A, Vol. 09, No. 24, pp. 2197-2206 (1994), https://doi.org/10.1142/S0217732394002057

[TC68]

“Exact Solution for an N-Molecule-Radiation-Field Hamiltonian”, M. Tavis and F. W. Cummings Phys. Rev. 170, 379 (1968) https://doi.org/10.1103/PhysRev.170.379

[SKFP16]

“TRIQS/CTHYB: A continuous-time quantum Monte Carlo hybridisation expansion solver for quantum impurity problems”, P. Seth, I. Krivenko, M. Ferrero and O. Parcollet, Comp. Phys. Comm. 200, March 2016, 274-284, http://dx.doi.org/10.1016/j.cpc.2015.10.023 (section 4.2)