{{ message }}

Instantly share code, notes, and snippets.

# marekyggdrasil/qm.py Secret

Last active Feb 26, 2022
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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