Code for https://mareknarozniak.com/2020/10/14/jordan-wigner-transformation/ tutorial article

qm.py

import numpy as np
from qutip import qeye, sigmax, sigmay, sigmaz, tensor, create, destroy

def Is(i, levels=2): return [qeye(levels) for j in range(0, i)]

def Sx(N, i): return tensor(Is(i) + [sigmax()] + Is(N - i - 1))
def Sy(N, i): return tensor(Is(i) + [sigmay()] + Is(N - i - 1))
def Sz(N, i): return tensor(Is(i) + [sigmaz()] + Is(N - i - 1))

def I(N): return Sz(N, 0)*Sz(N, 0)

def b_(N, Np, i): return tensor(Is(i, levels=Np+1) + [destroy(Np+1)] + Is(N - i - 1, levels=Np+1))
def bd(N, Np, i): return tensor(Is(i, levels=Np+1) + [create(Np+1)] + Is(N - i - 1, levels=Np+1))

def osum(lst): return np.sum(np.array(lst, dtype=object))

def oprd(lst, d=None):
    if len(lst) == 0:
        return d
    p = lst[0]
    for U in lst[1:]:
        p = p*U
    return p

# N is the exponent, L is the length of the chain
def opow(L, op, N): return oprd([op for i in range(N)])


def commutator(A, B):
    return A*B - B*A


def anticommutator(A, B):
    return A*B + B*A


def a_(N, n, Opers=None):
    Sa, Sb, Sc = Sz, Sx, Sy
    if Opers is not None:
        Sa, Sb, Sc = Opers
    return oprd([Sa(N, j) for j in range(n)], d=I(N))*(Sb(N, n) + 1j*Sc(N, n))/2.


def ad(N, n, Opers=None):
    Sa, Sb, Sc = Sz, Sx, Sy
    if Opers is not None:
        Sa, Sb, Sc = Opers
    return oprd([Sa(N, j) for j in range(n)], d=I(N))*(Sb(N, n) - 1j*Sc(N, n))/2.

test_operators.py

import pytest
import math
import itertools

from qm import Sx, Sy, Sz
from qm import a_, ad, b_, bd, I
from qm import opow
from qm import commutator, anticommutator


def fermionTestName(param):
    N, jw = param
    (_, la), (_, lb), (_, lc) = jw
    return 'N={0},JW={1}'.format(str(N), la + lb + lc)


Ns = [2, 3, 4]
Nps = [2, 3, 4]
jws = itertools.permutations([(Sx, 'X'), (Sy, 'Y'), (Sz, 'Z')])

fermion_params = list(itertools.product(Ns, jws))
fermion_names = [fermionTestName(param) for param in fermion_params]

boson_params = list(itertools.product(Ns, Nps))
boson_names = ['N={0},Ns={1}'.format(str(N), str(Ns)) for N, Ns in boson_params]


@pytest.mark.parametrize('N,jw', fermion_params, ids=fermion_names)
def testFermions(N, jw):
    (Sa, _), (Sb, _), (Sc, _) = jw
    Opers = Sa, Sb, Sc
    zero = 0.*I(N)
    # test all the pairs
    for n in range(N):
        a_n = a_(N, n, Opers=Opers)
        adn = ad(N, n, Opers=Opers)
        for np in range(N):
            a_np = a_(N, np, Opers=Opers)
            adnp = ad(N, np, Opers=Opers)
            assert anticommutator(a_n, a_np) == zero
            if n == np:
                assert anticommutator(a_n, adnp) == I(N)
            else:
                assert anticommutator(a_n, adnp) == zero
        assert a_n*a_n == zero
        assert adn*adn == zero


@pytest.mark.parametrize('N,Np', boson_params, ids=boson_names)
def testBosons(N, Np):
    zero = 0.*b_(N, Np, 0)*bd(N, Np, 0)
    # test all the pairs
    for n in range(N):
        b_n = b_(N, Np, n)
        bdn = bd(N, Np, n)
        for np in range(N):
            b_np = b_(N, Np, np)
            bdnp = bd(N, Np, np)
            # test anticommutation properties
            assert commutator(b_n, b_np) == zero
            LHS = commutator(b_n, bdnp)
            RHS = zero
            if n == np:
                NpF = math.factorial(Np)
                RHS = (1. - ((Np+1)/NpF)*opow(N, bdn, Np)*opow(N, b_n, Np))
            assert LHS == RHS