Introduction
HQS Quantum Solver enables you to run precise simulations of quantum systems by accounting for all interactions between the particles in a system. HQS Quantum Solver gives you easy access to high-performance routines and solvers needed to work with quantum systems.
Providing a consistent and extensible interface, it offers routines for the construction of Hamiltonians and operators in either full configurational space or subspaces characterized by conserved quantities such as the spin or the particle number.
Applications
HQS Quantum Solver is a powerful tool for researchers and enthusiasts in the field of condensed matter physics and beyond. It can be used for a wide range of applications, including but not limited to:
-
Simulating time-evolution of spin and fermionic systems.
-
Evaluating spectral functions of interacting fermionic systems and dynamic correlation function for fermionic and spin-spin problems.
-
Exploring quench scenarios in various lattice structures.
-
Studying interacting band structures for spinless and spinful fermion lattices using CPT.
Getting started
HQS Quantum Solver is a tool that can be used in conjunction with other HQStage
modules by HQS Quantum Simulations GmbH. To install this module run
hqstage install hqs-quantum-solver
To use HQS Quantum Solver, you furthermore need to install the Intel Math Kernel Library (MKL), which is described in the Installation section.
For a collection of examples to start using HQS Quantum Solver please refer to the Getting Started section.
Features
-
User-Friendly Interface: The HQS Quantum Solver provides an a user-friendly interface for creating lattice models. It also connects to the
HQS Spin Mapper
andstruqture
software packages by HQS Quantum Simulations GmbH, which allows for arbitrary user-defined models. Making it easy for users to construct Hamiltonians and operators effortlessly. -
Enables exploitation of symmetries like particle number, spin conservation, or just the fermion parity for BCS type models for superconductivity.
-
Extensibility: HQS Quantum Solver allows users to integrate their own backends, tools, and algorithms.
-
Interoperable: Integrates well with the scientific Python ecosystem, especially with NumPy and SciPy.
-
Calculation of static and dynamical correlations functions in frequency domain including the correction vector approach and the expansion in Chebyshev polynomials.
The API documentation can be found here.
Installation
HQS Quantum Solver
Visit cloud.quantumsimulations.de/software, search for "HQS Quantum Solver" in the package list and select "Download". Follow the instructions on that page to install HQS Quantum Solver.
To use HQS Quantum Solver, you furthermore need to install the Intel Math Kernel Library (MKL), which is discussed in the next section.
MKL
Installing MKL such that HQS Quantum Solver can access it can be done in three different ways.
-
You can use HQStage to install the MKL into the currently active virtual environment.
hqstage install mkl
-
You can manually provide a version of the MKL by making sure that the file
libmkl_rt.so
is found by the dynamic linker. That means, that a system-wide installation should be found automatically. -
You can install the MKL via pip. However, you need to create a symlink for the file
libmkl_rt.so
. The commands below perform the necessary steps.pip install mkl PLATLIB_DIR="$(python -c "import sysconfig; print(sysconfig.get_path('platlib'))")" ln -fs libmkl_rt.so.2 "${PLATLIB_DIR}/../../libmkl_rt.so"
Getting Started
Since quantum mechanical systems are described by operators acting on state vectors living in Hilbert spaces, this library is mainly concerned with providing such operators and functions that manipulate state vectors. In this chapter, you will learn the essentials of using the HQS Quantum Solver software by considering an example computation.
In this and the following examples, we shall consider fermion systems, but other types of systems are available as well.
Let us consider the tight-binding Hamiltonian defined on a chain lattice, i.e., the Hamiltonian is given by where is the number of lattice sites, the strength of the hopping term and the annihilation operator at site .

Imports and Parameters
First, we need to import the functions, classes and modules that are required for the script. If you want to follow along, copy the following imports into your script file. We will explain the classes and functions as we go along.
# Title : HQS Quantum Solver Getting Started Guide
# Filename : getting_started.py
# Python, NumPy, SciPy, and Matplotlib imports
from math import pi, sqrt
import numpy as np
from scipy.linalg import toeplitz, norm
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt
# HQS Quantum Solver imports
from hqs_quantum_solver.spinless_fermions import (
VectorSpace, Operator, hopping, number, annihilation
)
from hqs_quantum_solver.util import unit_vector
Next, we define the parameters that we will refer to in the script.
M = 16 # site count
N = 1 # particle count
IE = 0 # number of energy level to plot
With that out of the way, we can start setting up the tight-binding Hamiltonian.
Setting up the Hamiltonian
We want to look at a system given by the tight-binding Hamiltonian. Thus, the first step is to construct said Hamiltonian. Quantum Solver provides a generic interface for constructing quantum mechanical operators, as well as several convenience functions.
The first step in constructing an operator is to define the vector space
the operator will be acting on. In quantum mechanics we consider vector
spaces where each basis vector corresponds to a state of the system. Hence, to
define a Quantum Solver vector space, we need to specify the states of the
system. For this particular example, we want the states that have
exactly N
particles distributed across M
sites. The corresponding
vector space is constructed by creating an instance of the
VectorSpace
in the following way.
v = VectorSpace(sites=M, particle_numbers=N)
Note that Quantum Solver is also able to create vector spaces for a varying number of particles (see spinless_fermions.VectorSpace, HQS Quantum Solver Reference). Fixing the number of particles, however, reduces the amount of memory and computation time Quantum Solver needs.
The second step, is to define the mathematical terms that define the operator.
In Quantum Solver, operators are built by combining multiple mathematical
terms. For the current example, we use the function called
hopping
,
which is part of the
spinless_fermions
module, and which implements the term
where is the, so called, hopping matrix.
If we define the matrix to be
then the term in the previous formula becomes the tight-binding
Hamiltonian defined at the beginning of this chapter.
This matrix can be easily constructed using the
toeplitz
function from SciPy in combination with the
unit_vector
utility function.
t = 1
h0 = toeplitz(-t * unit_vector(v.sites, 1))
With this matrix at hand, we can construct the Hamiltonian, by passing
the term and the vector space defining the operator to the constructor of the
Operator
class.
H = Operator(hopping(h0), domain=v)
The Quantum Solver operators primarily provide matrix-vector and
matrix-matrix routines and are also directly compatible with most sparse
matrix routines from SciPy. More precisely, most operators, like the one we
have constructed here, implement the methods defined in the
OperatorProtocol
class, as well as the SciPy LinearOperator
interface. As a first example for using the
operator, we want to look at its spectrum. For that, we can directly use
the eigsh
function from SciPy.
k = min(12, H.shape[0] - 1) # Number of eigenvalues to compute
evals, evecs = eigsh(H, k=k, which="SA")
print(evals)
plt.figure("eigenvalues", clear=True)
plt.title(f"Eigenvalues (M = {M}, N = {N})")
plt.plot(evals, ".")
plt.xlabel("Number of Eigenvalue")
plt.ylabel("Eigenvalue")
plt.show()
Now, we have enough code to run the script for the first time. Doing so should produce the plot shown below, which shows the ten smallest eigenvalues of the Hamiltonian, and hence, the lowest energy levels of the system.

Before continuing, modify the values of and , run the script, and look at the result a few times.
Note, for a fixed number of particles N
the number of states
Quantum Solver has to store is
. Hence, the number of states can grow rather quickly when
increasing N
while M
is large. This problem gets worse when allowing
the particle number to vary. Keep in mind, while
Quantum Solver provides you with highly accurate computations, you will
have to provide enough computational resources for it to shine.
Below, we have listed the 5 smallest eigenvalues for M = 16
and different
values of N
.
# | |||
---|---|---|---|
1 | -1.97 | -3.83 | -5.53 |
2 | -1.86 | -3.67 | -5.31 |
3 | -1.70 | -3.57 | -5.14 |
4 | -1.48 | -3.44 | -5.04 |
5 | -1.21 | -3.34 | -5.04 |
We can make our first observation: Adding the first two values of the column gives you the first value in the column, and adding the first three values of the column will give you the first value in the column. The reason for this is that in a two- and three-particle system of non-interacting paricles, the lowest energy to obtain is the one of the particles occupying the lowest two and three energy levels, respectively. Since we are dealing with a fermionic system, the particles are not allowed to occupy the same state, which explains this behavior.
We shall later come back to this observation.
Note how only few lines of code were necessary to produce physically meaningful results.
Site Occupations
To take this example a little bit further, let us look at the computation of site occupation expectation values. More precisely, we want to compute for , where and
To evaluate this term, we first have to create the number operator for
each site. The
number
function, is
provided for this purpose. This function takes the named argument site
and returns the corresponding term.
The term can then be passed to the constructor of the
Operator
class,
as we have done for the construction of the Hamiltonian in the previous
section.
Once we have constructed the operator, we can use the
dot
method, to get .
Then we use the
numpy.vdot
function, to compute the scalar product. With this information, we can
write the code that measures the occupation expectation value.
Below, we define a function site_occupation
which computes the expectation
value for a given state vector and site. We then plot the expectation
value for every site and for the state vector corresponding to the
energy level IE
.
def site_occupation(psi, j):
number_operator = Operator(number(site=j), domain=v)
return np.vdot(psi, number_operator.dot(psi))
occupations = np.array([site_occupation(evecs[:, IE], j) for j in range(M)])
plt.figure("site_occupation", clear=True)
plt.title(f"Site Occupation (M = {M}, N = {N}, IE = {IE})")
plt.plot(occupations, "x-")
plt.axis((0, M - 1, 0, 1))
plt.xlabel("j")
plt.ylabel("Expectation Value")
plt.show()

The site occupations for M=16
, N=1
, IE=0
.

The site occupations for M=16
, N=2
, IE=0
.
Energy occupations
While looking at site occupations is certainly interesting, for most applications it is more useful to look at energy occupations instead.
Consider the matrix given by The columns of this matrix are the eigenvectors of the matrix . It is a well-known fact that if we define we can write the Hamiltonian as where are the energy levels of the system. Consequently, we can measure the expectation value of the occupation of the -th energy level, by computing Note that Hence, we just need to compute and then compute the square norm of this vector.
First, we have to construct the annihilation operator , which is
a linear combination of the annihilation operators . The
annihilation operator reduces the number of particles in the system by one.
In case an operator changes the particle number, in addition to specifying
the domain, the codomain needs to be given as well.
To construct the codomain vector space, we can use the
copy
method on the existing vector space to obtain a modified copy.
v_minus = v.copy(particle_number_change=-1)
The term that describes the annihilation operator is given by the
annihilation
function.
This function takes the keyword argument site
to specify the site at which
to apply the annihilation operator, but, conveniently, also
takes the coef
argument, which allows specifying a linear combination
of annihilation operators. Setting coef
to the desired linear combination
and then passing this term, the domain, and the codomain to the
constructor of the Operator class constructs the desired operator
. Once the operator is constructed, the computation of the
energy occupation is straightforward.
J, K = np.meshgrid(np.arange(M), np.arange(M), indexing='ij')
Q = sqrt(2) / sqrt(M+1) * np.sin((J+1) * (K+1) * pi / (M+1))
def energy_occupation(psi, k):
annihilation_operator = Operator(
annihilation(coef=Q[:,k].conj()), domain=v, codomain=v_minus
)
return norm(annihilation_operator.dot(psi))**2
ks = np.arange(M)
expectation = np.array([energy_occupation(evecs[:, IE], k) for k in ks])
plt.figure("energy", clear=True)
plt.title(f"Energy Occupation (M = {M}, N = {N}, IE = {IE})")
plt.plot(ks, expectation, "x-")
plt.xlabel("p")
plt.ylabel("Expectation Value")
plt.show()
Below, you can find the plots that this script produces for different parameter configurations.

The energy occupations for M=16
, N=2
, IE=0
.

The energy occupations for M=16
, N=2
, IE=1
.

The energy occupations for M=16
, N=2
, IE=2
.
In these figures, you can see how the fermions arrange in the different one-particle energy states to achieve the energy levels we have observed in the table in the Setting up the Hamiltonian section.
Note that if you do not know the eigenvectors of the matrix , you can
easily obtain them by using the
eigh
function of NumPy.
Summary
To conclude this introduction, let us sum up the main takeaways of this chapter.
- Quantum Solver provides classes to work with quantum mechanical operators.
- Operators are constructed by specifying a term and a vector space.
- For operators that change the particle number, the codomain must be given as well.
- Operators provide operations on vectors and matrices (as described in
the
OperatorProtocol
class), like thedot
method, which applies the operator to a state vector.
Complete Code
# Title : HQS Quantum Solver Getting Started Guide
# Filename : getting_started.py
# Python, NumPy, SciPy, and Matplotlib imports
from math import pi, sqrt
import numpy as np
from scipy.linalg import toeplitz, norm
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt
# HQS Quantum Solver imports
from hqs_quantum_solver.spinless_fermions import (
VectorSpace, Operator, hopping, number, annihilation
)
from hqs_quantum_solver.util import unit_vector
# ===== Parameters =====
M = 16 # site count
N = 1 # particle count
IE = 0 # number of energy level to plot
# ===== Creation of the Hamiltonian =====
v = VectorSpace(sites=M, particle_numbers=N)
t = 1
h0 = toeplitz(-t * unit_vector(v.sites, 1))
H = Operator(hopping(h0), domain=v)
# ===== Computing the Eigenvalues and Eigenvectors =====
k = min(12, H.shape[0] - 1) # Number of eigenvalues to compute
evals, evecs = eigsh(H, k=k, which="SA")
print(evals)
plt.figure("eigenvalues", clear=True)
plt.title(f"Eigenvalues (M = {M}, N = {N})")
plt.plot(evals, ".")
plt.xlabel("Number of Eigenvalue")
plt.ylabel("Eigenvalue")
plt.show()
# ===== Plotting the Site Occupations =====
def site_occupation(psi, j):
number_operator = Operator(number(site=j), domain=v)
return np.vdot(psi, number_operator.dot(psi))
occupations = np.array([site_occupation(evecs[:, IE], j) for j in range(M)])
plt.figure("site_occupation", clear=True)
plt.title(f"Site Occupation (M = {M}, N = {N}, IE = {IE})")
plt.plot(occupations, "x-")
plt.axis((0, M - 1, 0, 1))
plt.xlabel("j")
plt.ylabel("Expectation Value")
plt.show()
# ===== Plotting the Energy Occupations =====
v_minus = v.copy(particle_number_change=-1)
J, K = np.meshgrid(np.arange(M), np.arange(M), indexing='ij')
Q = sqrt(2) / sqrt(M+1) * np.sin((J+1) * (K+1) * pi / (M+1))
def energy_occupation(psi, k):
annihilation_operator = Operator(
annihilation(coef=Q[:,k].conj()), domain=v, codomain=v_minus
)
return norm(annihilation_operator.dot(psi))**2
ks = np.arange(M)
expectation = np.array([energy_occupation(evecs[:, IE], k) for k in ks])
plt.figure("energy", clear=True)
plt.title(f"Energy Occupation (M = {M}, N = {N}, IE = {IE})")
plt.plot(ks, expectation, "x-")
plt.xlabel("p")
plt.ylabel("Expectation Value")
plt.show()
Building Operators
Quantum Solver allows quantum mechanical operators to be constructed in
a flexible and convenient way. As we have seen in the
"Getting Started" section, quantum mechanical operators
are constructed by providing a mathematical term and a vector space.
The constructor of the
Operator
class,
however, does not only accept single terms, but arbitrary linear
combinations of terms, which makes this interface flexible.
To demonstrate this feature, let us consider the following example. Assume we want to construct the Hamiltonian , where . We start by performing the necessary imports and defining the Parameters.
# NumPy, SciPy, and Matplotlib imports
import numpy as np
from scipy.linalg import toeplitz
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt
# HQS Quantum Solver imports
from hqs_quantum_solver.spinless_fermions import (
VectorSpace,
Operator,
hopping,
interaction,
number,
)
from hqs_quantum_solver.util import unit_vector
# ===== Parameters =====
M = 16 # site count
N = 2 # particle count
t = 1 # hopping strength
U = 10 # interaction strength
Then, we create the
VectorSpace
object,
like in the "Getting Started" section, i.e.:
v = VectorSpace(sites=M, particle_numbers=N)
We can now consider the individual parts of the formula for the Hamiltonian.
We use the two functions
hopping
and
interaction
, which both
expect a matrix as input that describes between which sites a particle
can move and which pairs of sites expose particles to interactive forces,
respectively. The term
can be obtained by passing the matrix
to the hopping
function, and this matrix can be constructed by:
neighbors = toeplitz(unit_vector(v.sites, 1))
The term
can be obtained by passing the matrix
to the interaction
function, and this matrix can be constructed by:
right_neighbors = toeplitz(np.zeros(v.sites), unit_vector(v.sites, 1))
We can now combine the individual terms, to construct the desired Hamiltonian.
H = Operator(-t * hopping(neighbors) + U * interaction(right_neighbors), domain=v)
Using the same code as discussed in the "Getting Started" section, we can compute the expectation value of the occupation for different values of .
# ===== Computing the Eigenvalues and Eigenvectors =====
k = min(12, H.shape[0] - 1) # Number of eigenvalues to compute
evals, evecs = eigsh(H, k=k, which="SA")
groundstate = evecs[:, 0]
# ===== Plotting the Site Occupations =====
def site_occupation(psi, j):
n = Operator(number(site=j), domain=v)
return np.vdot(psi, n.dot(psi))
occupations = np.array([site_occupation(groundstate, j) for j in range(M)])
plt.figure("site_occupation", clear=True)
plt.title(f"Site Occupation (M = {M}, N = {N}, t={t}, U={U})")
plt.plot(occupations, "x-")
plt.axis((0, M - 1, 0, 1))
plt.xlabel("j")
plt.ylabel("Expectation Value")
plt.show()

The site occupations for .

The site occupations for .

The site occupations for .
Complete Code
# Title : HQS Quantum Solver Hubbard Interaction
# Filename : hubbard_interaction.py
# NumPy, SciPy, and Matplotlib imports
import numpy as np
from scipy.linalg import toeplitz
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt
# HQS Quantum Solver imports
from hqs_quantum_solver.spinless_fermions import (
VectorSpace,
Operator,
hopping,
interaction,
number,
)
from hqs_quantum_solver.util import unit_vector
# ===== Parameters =====
M = 16 # site count
N = 2 # particle count
t = 1 # hopping strength
U = 10 # interaction strength
# ===== Creation of the Hamiltonian =====
v = VectorSpace(sites=M, particle_numbers=N)
neighbors = toeplitz(unit_vector(v.sites, 1))
right_neighbors = toeplitz(np.zeros(v.sites), unit_vector(v.sites, 1))
H = Operator(-t * hopping(neighbors) + U * interaction(right_neighbors), domain=v)
# ===== Computing the Eigenvalues and Eigenvectors =====
k = min(12, H.shape[0] - 1) # Number of eigenvalues to compute
evals, evecs = eigsh(H, k=k, which="SA")
groundstate = evecs[:, 0]
# ===== Plotting the Site Occupations =====
def site_occupation(psi, j):
n = Operator(number(site=j), domain=v)
return np.vdot(psi, n.dot(psi))
occupations = np.array([site_occupation(groundstate, j) for j in range(M)])
plt.figure("site_occupation", clear=True)
plt.title(f"Site Occupation (M = {M}, N = {N}, t={t}, U={U})")
plt.plot(occupations, "x-")
plt.axis((0, M - 1, 0, 1))
plt.xlabel("j")
plt.ylabel("Expectation Value")
plt.show()
Lattice Builder
When building Hamiltonians with Quantum Solver, it is often required to create
adjacency matrices that represent the lattice structure of the system under
consideration, e.g., for the
hopping
term.
Building those terms in 1d is relatively straight forward. For higher
dimensions and more complicated lattices, however, this becomes more
difficult. The HQS Lattice Builder
and Lattice Validator
are libraries that simplify this process and can directly be used
with Quantum Solver.
As an example, we consider the computation of the groundstate of the tight-binding Hamiltonian, but this time on a rectangular lattice. The Hamiltonian is given by where is the set of all tuples such that is a neighbor of .

We start the example by specifying the required imports.
# NumPy, SciPy, and Matplotlib imports
import numpy as np
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt
# HQS Quantum Solver imports
from hqs_quantum_solver.spinless_fermions import (
VectorSpace, Operator, lattice_terms, annihilation
)
# HQS Lattice Builder/Validator imports
import lattice_builder
import lattice_validator
Then, we define the parameters of the simulation.
Mx = 5 # sites in x-direction
My = 4 # sitex in y-direction
N = 1 # particle count
t = 1 # hopping strengh
The next step is to create a Lattice Builder/Validator configuration
that specifies the Hamiltonian on an Mx
My
rectangular lattice.
In this example we have just one type of atom (lattice site) that is
repeated in the direction of the vectors in "lattice_vectors"
.
The hopping term is defined in the bonds
variable, which adds an entry
of value -t
into the hopping matrix for neighboring lattice sites.
For more details on the Lattice Builder configuration, take a look at the Hamiltonians in "Physics", HQS Lattice Documentation documentation, and at "Schema definitions", Lattice Validator Documentation.
system = {
"site_type": "spinless_fermions",
"cluster_size": [Mx, My, 1],
"system_boundary_conditions": ["hard-wall", "hard-wall", "hard-wall"]
}
atoms = [{"id": 0}]
bonds = [
{"id_from": 0, "id_to": 0, "translation": [1, 0, 0], "t": -t},
{"id_from": 0, "id_to": 0, "translation": [0, 1, 0], "t": -t},
]
unitcell = {
"atoms": atoms,
"bonds": bonds,
"lattice_vectors": [[1, 0, 0], [0, 1, 0]]
}
conf = {
"system": system,
"unitcell": unitcell
}
The configuration needs to be validated and normalized, which is
done by using the Lattice Validator. Then, a Lattice Builder instance
can be created. This instance is then passed to the
lattice_terms
function, to convert it into an operator definition.
Note, if you do not need access to the Lattice Builder instance, you can also
directly pass the conf
dict to lattice_terms
, which will take care of
the validation and builder construction.
is_valid, conf = lattice_validator.Validator().validateConfiguration(conf)
if not is_valid:
raise ValueError(f"Configuration not valid: {conf}")
builder = lattice_builder.Builder(conf)
v = VectorSpace(sites=Mx * My, particle_numbers=N)
H = Operator(lattice_terms(builder), domain=v)
After having obtained the Hamiltonian operator, the computation of the
groundstate and the definition of the site_occupation
function are
the same as in the "Getting Started" section.
k = min(12, H.shape[0] - 1) # Number of eigenvalues to compute
evals, evecs = eigsh(H, k=k, which="SA")
groundstate = evecs[:, 0]
v_minus = v.copy(particle_number_change=-1)
def site_occupation(psi, j):
annihilation_operator = Operator(annihilation(site=j), domain=v, codomain=v_minus)
c_psi = annihilation_operator.dot(psi)
return c_psi.conj() @ c_psi
To plot the result we use the get_atom_positions
method from the Lattice
Builder to show the groundstate occupation on a 2d grid.
occupations = np.zeros((My, Mx))
for j, position in enumerate(builder.get_atom_positions().tolist()):
x, y = round(position[0]), round(position[1])
# Set value for row=y and column=x.
occupations[y, x] = site_occupation(groundstate, j)
plt.figure("site_occupation", clear=True)
plt.title(f"Site Occupation (lattice: {Mx}x{My}, particles: {N})")
plt.imshow(occupations, origin="lower", aspect="equal")
plt.xlabel("x")
plt.ylabel("y")
plt.xticks(np.arange(Mx))
plt.yticks(np.arange(My))
plt.colorbar()
plt.show()
Finally, we can run the code, to obtain a plot of the occupation, as shown below.

Complete Code
# Title : HQS Quantum Solver using the Lattice Builder
# Filename : getting_started.py
# NumPy, SciPy, and Matplotlib imports
import numpy as np
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt
# HQS Quantum Solver imports
from hqs_quantum_solver.spinless_fermions import (
VectorSpace, Operator, lattice_terms, annihilation
)
# HQS Lattice Builder/Validator imports
import lattice_builder
import lattice_validator
# ===== Parameters =====
Mx = 5 # sites in x-direction
My = 4 # sitex in y-direction
N = 1 # particle count
t = 1 # hopping strengh
# ===== Definition of the Lattice =====
system = {
"site_type": "spinless_fermions",
"cluster_size": [Mx, My, 1],
"system_boundary_conditions": ["hard-wall", "hard-wall", "hard-wall"]
}
atoms = [{"id": 0}]
bonds = [
{"id_from": 0, "id_to": 0, "translation": [1, 0, 0], "t": -t},
{"id_from": 0, "id_to": 0, "translation": [0, 1, 0], "t": -t},
]
unitcell = {
"atoms": atoms,
"bonds": bonds,
"lattice_vectors": [[1, 0, 0], [0, 1, 0]]
}
conf = {
"system": system,
"unitcell": unitcell
}
# ===== Creation of the Hamiltonian =====
is_valid, conf = lattice_validator.Validator().validateConfiguration(conf)
if not is_valid:
raise ValueError(f"Configuration not valid: {conf}")
builder = lattice_builder.Builder(conf)
v = VectorSpace(sites=Mx * My, particle_numbers=N)
H = Operator(lattice_terms(builder), domain=v)
# ===== Computing the Eigenvalues and Eigenvectors =====
k = min(12, H.shape[0] - 1) # Number of eigenvalues to compute
evals, evecs = eigsh(H, k=k, which="SA")
groundstate = evecs[:, 0]
# ===== Plotting the Site Occupations =====
v_minus = v.copy(particle_number_change=-1)
def site_occupation(psi, j):
annihilation_operator = Operator(annihilation(site=j), domain=v, codomain=v_minus)
c_psi = annihilation_operator.dot(psi)
return c_psi.conj() @ c_psi
occupations = np.zeros((My, Mx))
for j, position in enumerate(builder.get_atom_positions().tolist()):
x, y = round(position[0]), round(position[1])
# Set value for row=y and column=x.
occupations[y, x] = site_occupation(groundstate, j)
plt.figure("site_occupation", clear=True)
plt.title(f"Site Occupation (lattice: {Mx}x{My}, particles: {N})")
plt.imshow(occupations, origin="lower", aspect="equal")
plt.xlabel("x")
plt.ylabel("y")
plt.xticks(np.arange(Mx))
plt.yticks(np.arange(My))
plt.colorbar()
plt.show()
System Types
Up until this point, we have only considered quantum mechanical systems
composed of "spinless fermions". HQS Quantum Solver, however, provides
further types of quantum systems. Each of these systems is implemented
in its own module, and all these modules are structured in a similar way,
meaning that they all provide a VectorSpace
, an Operator
, and a set
of functions providing operator terms.
Spins
The spins
module implements systems
consisting of spin- particles.
The VectorSpace
is spanned by
vectors of the form
As a shorthand, we usually write, e.g.,
With this module, e.g., the
quantum Heisenberg model,
can be simulated. See the spin_waves.ipynb
example for details.
Spinless Fermions
The spinless_fermions
module implements
systems of fermions having a site index but no (dedicated) spin index.
The VectorSpace
is spanned
by vectors of the form
and is the creation operator for site .
With this module, e.g., the tight-binding model, can be simulated, as was shown in the "Getting Started" section.
Spinful Fermions
The spinful_fermions
module implements
system of fermions having a site index and a spin index.
The VectorSpace
is spanned
by vectors of the form
and is the creation operator for site and spin
polarization .
With this module, e.g., the Hubbard model,
can be simulated. See the spinful_fermions_groundstate.ipynb
example for
details.
Bosons
The bosons
module implements systems
of bosons.
The VectorSpace
is spanned
by vectors of the form
and is the creation operator for site .
The variable is the maximal occupation of site .
With this module, e.g., the Bose-Hubbard model,
can be simulated. See the bose_hubbard_groundstate.ipynb
example for
details.
Examples
Before looking at the examples, it is advised that you read the Getting Started Guide.
The examples described below can be accessed via HQStage.
-
spinful_fermions_groundstate.ipynb
In this example, we show how to compute the ground state of the Hamiltonian of a many-particle system.
-
spin_waves.ipynb
This example introduces spin systems and time evolution. We consider the time evolution of a Heisenberg spin model in an excited state.
-
nonequilibrium_electron_transport.ipynb
This example simulates the time evolution of a system of electrons starting in a nonequilibrium state.
-
spinless_fermions_greens_function.ipynb
Example for the Calculation of the Spectral function using HQS Quantum Solver and Lattice Functions.
This method combines operators, Chebyshev expansion, cluster perturbation theory and clever Green's function mathematics to create the reciprocal-space resolved spectral function (a. k. a. the interacting band structure) in a honeycomb lattice causing an opening of the Dirac point.
-
spinful_fermions_greens_function.ipynb
Prototype example for the Calculation of the Spectral function of the Hubbard model using HQS Quantum Solver. In this approach we represent the wavefunction as a dyadic product of up and down states.
Similar to the previous example, this example computes the interacting band structure of a square lattice under the influence of interaction. We retrieve the expected behaviour of opening a gap at half-filling.
-
bose_hubbard_groundstate.ipynb
In this example we compute the groundstates of Bose-Hubbard systems, which is a simple bosonic system.
-
lidar.ipynb
In this scientific example we study properties of light emitted by multiple atoms, a fundamental model in quantum optics. We use the HQS Quantum Solver to time evolve the system according to the Lindblad equation of motion and calculate various expectation values of system operators describing properties of the emitted light.
API Reference
You can find the API reference here.
Custom Operators
Many algorithms in HQS Quantum Solver can be used with your own operator implementations. All that is required is that your implementation fulfills the requirements of the protocols used by the algorithm of interest.
Protocols
Protocols are a Python language feature that allows describing the requirements an object has to fulfill for being used in a particular part of the program. Whenever you encounter a protocol class in the type hints of a function or method, the protocol class defines which methods and properties the corresponding argument object is expected to have and the corresponding return type it has to provide, respectively. It is not necessary that your class inherits from the protocol class.
For example, the
evolve_krylov
function
expects a operator
argument which adheres to the
SimpleOperatorProtocol
.
To be allowed to pass an object as operator
, the object just needs
to implement all the methods and properties listed in
SimpleOperatorProtocol
,
which are the shape
and dtype
attributes and the dot
method.
Example
In this example we want to simulate the time evolution of a single spin in a magnetic field for different field stengths, i.e, we consider the Hamiltonian for different values of and measure the expectation value of . Instead of building a new operator for every value of we are interested in, we create a Hamiltonian class, that applies the correct value of on the fly. Avoiding recreation of operators can sometimes reduce the time of the computation. For a single-spin system, however, we will not see any benefit, because the system is too small. Nevertheless, we stick with this simple example for demonstration purposes.
We start by performing the necessary imports.
# Title : HQS Quantum Solver Larmor Precession
# Filename : larmor_precession.py
import numpy as np
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt
from hqs_quantum_solver.spins import VectorSpace, Operator, spin_x, spin_z
from hqs_quantum_solver.evolution import evolve_krylov, observers
Next, we construct the and operators, and we set the initial state to an eigenvalue of the operator.
v = VectorSpace(sites=1, total_spin_z="all")
Sx = Operator(spin_x(site=0), domain=v)
Sz = Operator(spin_z(site=0), domain=v)
eigvals, eigvecs = eigsh(Sx, k=1)
state = eigvecs[:, 0]
We can now define the Hamiltonian
class as a custom operator.
In this class, we need to define the shape
and dtype
property and
the dot
method. Note that the implementation of the dot
method
is required to have a parameter called out
that when not equal to
None
needs to be filled with the result of the operator application.
class Hamiltonian:
def __init__(self, B):
self.B = B
@property
def shape(self):
return Sz.shape
@property
def dtype(self):
return Sz.dtype
def dot(self, x, out=None):
result = -self.B * Sz.dot(x)
if out is not None:
np.copyto(out, result)
return result
We can now pass instances of the Hamiltonian class to the
evolve_krylov
function, like
any other operator provided by HQS Quantum Solver.
times = np.linspace(0, 1, 100)
Bs = [1, 2, 5, 10]
observations = np.empty((len(Bs), times.size))
for n in range(len(Bs)):
B = Bs[n]
result = evolve_krylov(Hamiltonian(B=B), state, times, observer=observers.expectation([Sx]))
observations[n, :] = np.array(result.observations)[:, 0]
plt.plot(
times, observations[0, :],
times, observations[1, :],
times, observations[2, :],
times, observations[3, :],
)
plt.legend([f"$B = {B}$" for B in Bs])
plt.xlabel(r"$t$")
plt.ylabel(r"$\langle S^x \rangle$")
plt.show()

The expectation value of measuring .
Complete Code
# Title : HQS Quantum Solver Larmor Precession
# Filename : larmor_precession.py
import numpy as np
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt
from hqs_quantum_solver.spins import VectorSpace, Operator, spin_x, spin_z
from hqs_quantum_solver.evolution import evolve_krylov, observers
v = VectorSpace(sites=1, total_spin_z="all")
Sx = Operator(spin_x(site=0), domain=v)
Sz = Operator(spin_z(site=0), domain=v)
eigvals, eigvecs = eigsh(Sx, k=1)
state = eigvecs[:, 0]
class Hamiltonian:
def __init__(self, B):
self.B = B
@property
def shape(self):
return Sz.shape
@property
def dtype(self):
return Sz.dtype
def dot(self, x, out=None):
result = -self.B * Sz.dot(x)
if out is not None:
np.copyto(out, result)
return result
times = np.linspace(0, 1, 100)
Bs = [1, 2, 5, 10]
observations = np.empty((len(Bs), times.size))
for n in range(len(Bs)):
B = Bs[n]
result = evolve_krylov(Hamiltonian(B=B), state, times, observer=observers.expectation([Sx]))
observations[n, :] = np.array(result.observations)[:, 0]
plt.plot(
times, observations[0, :],
times, observations[1, :],
times, observations[2, :],
times, observations[3, :],
)
plt.legend([f"$B = {B}$" for B in Bs])
plt.xlabel(r"$t$")
plt.ylabel(r"$\langle S^x \rangle$")
plt.show()