Skip to content

Instantly share code, notes, and snippets.

@marekyggdrasil

marekyggdrasil/qm.py Secret

Last active Oct 14, 2020
Embed
What would you like to do?
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 osum([Sz(N, i)*Sz(N, i) for i in range(N)])/N
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): return np.prod(np.array(lst, dtype=object))
def opow(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)])*(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)])*(Sb(N, n) - 1j*Sc(N, n))/2.
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(bdn, Np)*opow(b_n, Np))
assert LHS == RHS
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.