Tuto 1: first steps with the package

Dr. Saad Yalouz - Laboratoire de Chimie Quantique de Strasbourg, France - July 2022

Philosophy of the package

The philosophy of the QuantNBody package is to facilitate the implementation and manipulation of quantum many-body systems composed of electrons or bosons. To achieve this goal, the package has been designed to provide a fast and easy way to construct many-body operators and wave functions in a given many-body basis. In this way, it becomes possible to access quantities/objects of interest in a few lines of python to speed up and facilitate methodology development.

To that purpose, the package works with two fundamental ingredients.

  1. The first is the creation of a many-body basis (based on a total number of quantum particles and modes/orbitals to be filled) in which each operator can be represented.

  2. The second is the creation of the set of \(a^\dagger_p a_q\) hopping operators that are needed to construct any many-body operator that conserves the number of particles.

Once these two ingredients are created, the user can use predefined functions to construct different types of many-body operators (e.g. Hamiltonians, spin operators), and manipulate/view many-body quantum states. Note that the QuantNBody package has also been designed to provide flexibility for users to also create their own operators and functions based on the provided tools.

Nota Bene: For sake of simplicity, we will focus in these tutorials on fermionic systems.

Let us first import the package !

import quantnbody as qnb   # <==== General import
import numpy as np

Building a many-body basis

To build a many-body basis for a fermionic system, the QuantNBody package generates a list of many-body states which describe the repartition of \(N_{elec}\) electrons in \(2N_{MO}\) spin-orbitals. These states are numerically referenced by a list of kappa indices such that

\[\Big\lbrace |\kappa \rangle \Big\rbrace_{\textstyle \kappa=1}^{\textstyle \dim_H}\]

The dimension \(\dim_H\) of the many-body basis depends on the number of electron \(N_{elec}\) and spatial orbital \(N_{MO}\) via a binomial law such that

\[\dim_H = \binom{2N_{MO}}{N_{elec}}\]

A little example with \(N_{MO}=N_{elec}=2\): in this case, we should have 6 many-body states.

N_MO = N_elec = 2 # We define the number of MO and electrons

nbody_basis = qnb.fermionic.tools.build_nbody_basis( N_MO, N_elec ) # Building the nbody_basis

print('Shape  of the kappa states')
for s in range(len(nbody_basis)):
    print('| kappa={} >'.format(s), '=', nbody_basis[s])

Shape  of the kappa states
| kappa=0 > = [1 1 0 0]
| kappa=1 > = [1 0 1 0]
| kappa=2 > = [1 0 0 1]
| kappa=3 > = [0 1 1 0]
| kappa=4 > = [0 1 0 1]
| kappa=5 > = [0 0 1 1]

What is the meaning of these six bit strings ?

Here, each bit string represents a many-body state. As an example, let us check the first state for which we have

\[| \kappa = 0\rangle = | \underbrace{ \overbrace{1}^{ \textstyle {\alpha}}, \; \; \;\overbrace{1}^{ \textstyle {\beta}},}_{\textstyle 1st \ MO}\; \; \underbrace{\overbrace{0}^{ \textstyle {\alpha}}, \; \; \; \overbrace{0}^{ \textstyle {\beta}}}_{\textstyle 2nd \ MO} \rangle\]

Here we choose to structure the occupation numbers as follows

  • Each couple of terms refer to a same spatial orbital

  • Even indices refer to \(\alpha\)-spin-orbitals

  • Odd indices refer to \(\beta\)-spin-orbitals

Note

Building a matrix representation of a many-body operator, what does it mean ? For each configuration, we associate a unique \(\kappa\) index which defines a unique “numerical” vector. In practice, any numerical representation of a given many-body operator will be given in the numerical many-body basis indexed by the \(\kappa\). As an example, let us imagine we want to encode numerically a second quantization operator \(O\). In practice, this means that we create a matrix representation of this operator in the many-body basis such that

\[O = \sum_{\kappa, \kappa' =1}^{\dim_H} \langle \kappa' | O | \kappa \rangle \; | \kappa' \rangle\langle \kappa |\]

In practice, this indexing is realized by the QuantNBody package and then used as a central tool to build every matrix element of a given many-body operators.

Building and storing the \(a^\dagger_{p,\sigma} a_{q,\tau}\) operators

Once the list of many-body states is created, the next crucial point in the QuantNBody package consists in building the \(a^\dagger_{p,\sigma} a_{q,\tau}\) many-body operators.

In practice, these operators play a central role in many cases of study as soon as we have to deal with systems that are particle-number conserving. In this case, one can show that many objects (i.e. excitation operators, spin operators, reduced density matrices …) are built in practice using series of \(a^\dagger_{p,\sigma} a_{q,\tau}\) operators.

With the QuantNBody package, we build the \(a^\dagger_{p,\sigma} a_{q,\tau}\) operators once and for all and store them via a very simple command line. In this way, we will be able to use them later on for any type of developments.

The command line is simple and only requires the list of many-body states we built previously :

a_dagger_a = qnb.fermionic.tools.build_operator_a_dagger_a( nbody_basis )

How to get access to these operators once stored ?

The way each operator is stored follows the way we order the spin-orbitals in our many-body states. As an illustrative example, taking the following elements will return the associated many-body operators :

a_dagger_a[0,0] \(\longrightarrow a^\dagger_{0,\alpha} a_{0,\alpha}\)

a_dagger_a[1,0] \(\longrightarrow a^\dagger_{0,\beta} a_{0,\alpha}\)

a_dagger_a[10,1] \(\longrightarrow a^\dagger_{5,\alpha} a_{0,\beta}\)

In practice, the resulting many-body operators we get access to are expressed in the numerical many-body basis (i.e. the \(\kappa\) indices). To better illustrate that, let us consider a_dagger_a[0,0]. As explained before, this numerical object encodes \(a^\dagger_{0,\alpha} a_{0,\alpha}\) which can theoretically be expressed in the many-body basis as follows

\[a^\dagger_{0,\alpha} a_{0,\alpha} = \sum_{\kappa, \kappa' =1}^{\dim_H} \langle \kappa' | a^\dagger_{0,\alpha} a_{0,\alpha} | \kappa \rangle \; | \kappa' \rangle\langle \kappa |,\]

As the effect of this operator is simply to count the number of electrons in the first spin-orbital, we can simplify the previous expression as follows

\[a^\dagger_{0,\alpha} a_{0,\alpha} = | 1100 \rangle\langle 1100 | + | 1010 \rangle\langle 1010 | + | 1001 \rangle\langle 1001 |.\]

In terms of encoding in the numerical \(\kappa\) many-body basis, this operator can be rewritten as follows

\[a^\dagger_{0,\alpha} a_{0,\alpha} = | \kappa=0 \rangle\langle \kappa=0 | + | \kappa=1 \rangle\langle \kappa=1 | + | \kappa=2 \rangle\langle \kappa=2 |.\]

With QuantNBody, the object a_dagger_a[0,0] can be seen as the table (i.e. matrix) encoding all the matrix elements \(\langle \kappa' | a^\dagger_{0,\alpha} a_{0,\alpha} | \kappa \rangle\) of the associated many-body operator \(a^\dagger_{0,\alpha} a_{0,\alpha}\). This encoding is realized in practice via a sparse matrix representation which can be made dense as shown below

print(  "Sparse representation of a_dagger_a[0,0]" )
print( a_dagger_a[0,0] )

print( )
print( "Dense representation of a_dagger_a[0,0]" )
print( a_dagger_a[0,0].A )

Sparse representation of a_dagger_a[0,0]
  (0, 0)    1.0
  (1, 1)    1.0
  (2, 2)    1.0

Dense representation of a_dagger_a[0,0]
[[1. 0. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]]

Here, we see here that the effect of this operator is simply to count the number of electrons in the first spin-orbital. This explains why we only have ones on the three first elements of the diagonal !

Building our first many-body Hamiltonian : a fermi-Hubbard molecule

In this final part of the tutorial we will use the previously built a_dagger_a variable to implement a fermi-Hubbard molecule. In the local site basis, the model Hamiltonian is usually expressed such that:

\[\hat{H} = \color{blue}{\sum_{\langle i,j \rangle}^{N_{MO}} -t_{ij} \sum_{\sigma=\uparrow,\downarrow} (\hat{a}^\dagger_{j,\sigma}\hat{a}_{i,\sigma}+\hat{a}^\dagger_{i,\sigma}\hat{a}_{j,\sigma})} + \color{red}{\sum_i^{N_{MO}} \mu_{ii} \sum_{\sigma=\uparrow,\downarrow} \hat{a}^\dagger_{i,\sigma}\hat{a}_{i,\sigma} } + \color{black}{ \sum_i^{N_{MO}} U_{iiii} \hat{a}^\dagger_{i,\uparrow}\hat{a}_{i,\uparrow} \hat{a}^\dagger_{i,\downarrow}\hat{a}_{i,\downarrow} }\]
with :

  • \(t_{ij}\) the hopping terms between the pair of connected sites \(\langle i, j \rangle\).

  • \(\mu_{ii}\) the local chemical potential on site “\(i\)”.

  • \(U_{iiii}\) the local coulombic repulsion on site “\(i\)”.

We illustrate the shape of the system below
_images/graph.png

In a more general basis (not necessarily local) we have

\[\hat{H} =\sum_{\langle p,q \rangle}^{N_{MO}} -h_{pq} \sum_{\sigma=\uparrow,\downarrow} (\hat{a}^\dagger_{p,\sigma}\hat{a}_{q,\sigma}+\hat{a}^\dagger_{q,\sigma}\hat{a}_{p,\sigma}) + \sum_i^{N_{MO}} U_{p,q,r,s} \hat{a}^\dagger_{p,\uparrow}\hat{a}_{q,\uparrow} \hat{a}^\dagger_{r,\downarrow}\hat{a}_{s,\downarrow} ,\]

where, for commodity, we have introduced the one-body integrals \(h_{pq}\) which embed the hopping terms and the chemical potentials such as

\[h_{pq} = \sum_{i,j}^{N_{MO}} (-t_{ij} + \delta_{ij}\mu_{ii}) C_{i,p} C_{j,q},\]

and the “delocalized version” of the coulombic repulsion term

\[U_{pqrs} = \sum_{i}^{N_{MO}} U_{i,i,i,i} C_{i,p} C_{i,q} C_{i,r} C_{i,s},\]

where the matrix \({\bf C}\) encodes the Molecular Orbital coefficients (used if we want for example to express the Hamiltonian in a delocalized basis).

Building the Hamiltonian : To initiate the construction of the matrix representation of the operator in the many-body basis, we first define the hopping term \(t\) between the sites, the chemical potentials \(\mu\) and the electronic repulsion \(U\).

# Setup for the simulation ========
N_MO   = N_elec = 2
t_  = np.zeros((N_MO,N_MO))
U_  = np.zeros((N_MO,N_MO,N_MO,N_MO))
Mu_ = np.zeros((N_MO,N_MO))
for i in range(N_MO):
    U_[i,i,i,i]  =  1 * (1+i)  # Local coulombic repulsion
    Mu_[i,i]     = -1 * (1+i)  # Local chemical potential

    for j in range(i+1,N_MO):
        t_[i,j] = t_[j,i] = - 1  # hopping

h_ = t_  + np.diag( np.diag(Mu_) ) # Global one-body matrix = hoppings + chemical potentials

print( 't_=\n',t_ ,'\n')

print( 'Mu_=\n',Mu_ ,'\n')

print( 'h_=\n',h_ ,'\n')

t_=
 [[ 0. -1.]
 [-1.  0.]]

Mu_=
 [[-1.  0.]
 [ 0. -2.]]

h_=
 [[-1. -1.]
 [-1. -2.]]

To build the Hamiltonian, we simply have to pass the three following ingredients to an already built function:

  • Parameters of the model

  • The Many-body basis

  • The \(a^\dagger a\) operators

as shown below

H_fermi_hubbard = qnb.fermionic.tools.build_hamiltonian_fermi_hubbard( h_,
                                                                       U_,
                                                                       nbody_basis,
                                                                       a_dagger_a )

Similarily to the \(a^\dagger a\) operators, the Hamiltonian \(H\) is represented in the many-body basis with a native sparse representation (which can be made dense):

print('H (SPARSE) =' )
print(H_fermi_hubbard)

print()
print('H (DENSE) =' )
print(H_fermi_hubbard.A)

H (SPARSE) =
  (0, 0)    -1.0
  (0, 2)    -1.0
  (0, 3)    1.0
  (1, 1)    -3.0
  (2, 0)    -1.0
  (2, 2)    -3.0
  (2, 5)    -1.0
  (3, 0)    1.0
  (3, 3)    -3.0
  (3, 5)    1.0
  (4, 4)    -3.0
  (5, 2)    -1.0
  (5, 3)    1.0
  (5, 5)    -2.0

H (DENSE) =
[[-1.  0. -1.  1.  0.  0.]
 [ 0. -3.  0.  0.  0.  0.]
 [-1.  0. -3.  0.  0. -1.]
 [ 1.  0.  0. -3.  0.  1.]
 [ 0.  0.  0.  0. -3.  0.]
 [ 0.  0. -1.  1.  0. -2.]]

Once \(H\) is built, we can diagonalize the resulting matrix using for example the numpy library.

eig_energies, eig_vectors =  np.linalg.eigh(H_fermi_hubbard.A)

print('Energies =', eig_energies[:4] )

Energies = [-4.41421356 -3.         -3.         -3.        ]

And finally, we can call a very useful function from the QuantNBody package that help visualizing the shape of a wavefunction as shown below. This function lists the most important many-body states contributing to the wavefunction with the associated coefficients in front.

WFT_to_analyse = eig_vectors[:,0]

# Visualizing the groundstate in the many-body basis
qnb.fermionic.tools.visualize_wft( WFT_to_analyse, nbody_basis ) # <=== FCT IN THE PACKAGE
print()

-----------
 Coeff.     N-body state and index
-------     ----------------------
-0.57454        |0110⟩        #3
+0.57454        |1001⟩        #2
+0.47596        |0011⟩        #5
+0.33656        |1100⟩        #0