dylanljones / exactdiag

Exact diagonalization of fermionic (many-body) systems

Exact diagonalization

⚠️	This project is still under heavy development and might contain bugs or have breaking API changes in the future.

🔧 Installation

Install via pip from github:

pip install git+https://github.com/dylanljones/exactdiag.git@VERSION

or download/clone the package, navigate to the root directory and install via

pip install .

🚀 Quick-Start

Basis

A Basis object can be initalized with the number of sites in the (many-body) system:

import exactdiag as ed

basis = ed.Basis(num_sites=3)

The corresponding states of a particle sector can be obtained by calling:

sector = basis.get_sector(n_up=1, n_dn=1)

If no filling for a spin-sector is passed all possible fillings are included. The labels of all states in a sector can be created by the state_labels method:

>>> sector.state_labels()
['..⇅', '.↓↑', '↓.↑', '.↑↓', '.⇅.', '↓↑.', '↑.↓', '↑↓.', '⇅..']

The states of a sector can be iterated by the states-property. Each state consists of an up- and down-SpinState:

state = list(sector.states)[0]
up_state = state.up
dn_state = state.dn

Each SpinState provides methods for optaining information about the state, for example:

>>> up_state.binstr(width=3)
001
>>> up_state.n
1
>>> up_state.occupations()
[1]
>>> up_state.occ(0)
1
>>> up_state.occ(1)
0
>>> up_state.occ(2)
0

Operators

The operators-module provides the base-class LinearOperator based on scipy.LinearOperator. A simple sparse implementation of a Hamiltonian is also included.

import exactdiag as ed

size = 5
rows = [0, 1, 2, 3, 4, 0, 1, 2, 3, 1, 2, 3, 4]
cols = [0, 1, 2, 3, 4, 1, 2, 3, 4, 0, 1, 2, 3]
data = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1]
indices = (rows, cols)
hamop = ed.HamiltonOperator(size, data, indices)

Converting the operator to an array yields

>>> hamop.array()
[[0 1 0 0 0]
 [1 0 1 0 0]
 [0 1 0 1 0]
 [0 0 1 0 1]
 [0 0 0 1 0]]

Many-Body Hamiltonian matrices can be constructed by projecting the elements onto a basis sector. First, the basis has to be initialized:

import exactdiag as ed

basis = ed.Basis(num_sites=2)
sector = basis.get_sector()  # Full basis

The Hubbard Hamilton-operator, for example, can then be constructed as follows:

def hubbard_hamiltonian_data(sector):
    up_states = sector.up_states
    dn_states = sector.dn_states
    # Hubbard interaction
    yield from ed.project_hubbard_inter(up_states, dn_states, u=[2.0, 2.0])
    # Hopping between sites 0 and 1
    yield from ed.project_hopping(up_states, dn_states, site1=0, site2=1, hop=1.0)

rows, cols, data = list(), list(), list()
for i, j, val in hubbard_hamiltonian_data(sector):
    rows.append(i)
    cols.append(j)
    data.append(val)

hamop = ed.HamiltonOperator(sector.size, data, (rows, cols))

Models

Some methods require a model object to work. Users can define their own models or use one of the following included models:

Module	Description	Lattice support
abc	Model-Parameter container and abstract base classes	-
anderson	Anderson impurity models	❌
hubbard	Hubbard model	✔️

A custom model can be defined by inheriting from the abstract base classes. Many-body models, for example, have to implement the _hamiltonian_data method, which generates the rows, columns and values of the Hamilton operator for each basis sector:

import exactdiag as ed


class CustomModel(ed.models.AbstractManyBodyModel):

    def __init__(self, num_sites, eps=0.0, ...):
        super().__init__(num_sites, eps=eps, ...)

    def _hamiltonian_data(self, up_states, dn_states):
        yield from ed.project_onsite_energy(up_states, dn_states, self.eps)
        ...

The Hamilton operator can then be acessed for the full basis or a specific sector in operator or matrix representation:

model = ed.models.HubbardModel(num_sites=2, neighbors=[[0, 1]], inter=2)

# Sector n↑=1, n↓=1
hamop = model.hamilton_operator(n_up=1, n_dn=1)

# Full basis
sector = model.basis.get_sector()
ham = model.hamiltonian(sector=sector)
ed.matshow(ham, ticklabels=sector.state_labels(), values=True)

Green's function

Using a custom defined model or one of the included models the one-particle (many-body) Green's function can be computed, for example using the Lehmann sum:

import numpy as np
import exactdiag as ed

model = ed.models.HubbardModel.chain(num_sites=7, inter=4.0, beta=10.0).hf()

z = np.linspace(-8, +8, 1001) + 1e-1j
gf = ed.gf_lehmann(model, z, i=3, sigma=ed.UP)[0]

exactdiag's Issues

Numpy/Scipy import error

Running the tests in a GitHub Workflow on Windows using Python 3.10 results in Numpy/Scipy import errors:

ImportError while importing test module 'D:\a\exactdiag\exactdiag\tests\test_basis.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
C:\hostedtoolcache\windows\Python\3.10.7\x64\lib\importlib\__init__.py:1[26](https://github.com/dylanljones/exactdiag/actions/runs/3226005826/jobs/5279006230#step:6:27): in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
tests\test_basis.py:11: in <module>
    from exactdiag import basis
exactdiag\__init__.py:[27](https://github.com/dylanljones/exactdiag/actions/runs/3226005826/jobs/5279006230#step:6:28): in <module>
    from .matrix import (
exactdiag\matrix.py:11: in <module>
    from scipy.sparse import spmatrix
C:\hostedtoolcache\windows\Python\3.10.7\x64\lib\site-packages\scipy\sparse\__init__.py:267: in <module>
    from ._csr import *
C:\hostedtoolcache\windows\Python\3.10.7\x64\lib\site-packages\scipy\sparse\_csr.py:10: in <module>
    from ._sparsetools import (csr_tocsc, csr_tobsr, csr_count_blocks,
E   ImportError: numpy.core.multiarray failed to import
------------------------------- Captured stderr -------------------------------
RuntimeError: module compiled against API version 0xf but this version of numpy is 0xe

Running the tests local with the same dependencies from the requirements.txt file works as expected for Python versions 3.7-3.10

Green's function using Lehmann sum does not return the expected result for U=0 and N>2

Describe the bug
The Lehmann Green's function does not match the non-interacting GF for U=0 and bulk sites (not 0 or N-1)

To Reproduce
Steps to reproduce the behavior:

import numpy as np
import matplotlib.pyplot as plt
import exactdiag as ed


def greens0(model, z):
    nsites = model.num_sites
    ham = np.zeros((nsites, nsites), dtype=np.float64)
    np.fill_diagonal(ham, model.eps)
    for i, j in model.neighbors:
        ham[i, j] = ham[j, i] = model.hop

    return ed.gf0_pole(ham, z=z, mode="full")


model = ed.models.HubbardModel.chain(3, inter=0.0, beta=200.0).hf()
z = np.linspace(-5, +5, 1001) + 1e-1j

gfmat = ed.greens.compute_gf_full(model, z)
gfmat0 = greens0(model, z)

fig, ax = plt.subplots()
ax.plot(z.real, -gfmat0[:, 1, 1].imag)
ax.plot(z.real, -gfmat[:, 1, 1].imag, ls="--")
plt.show()

Expected behavior
The Green's function should match the non-interacting GF