%matplotlib inline
import numpy as np
import scipy.sparse.linalg as la
Cout how many bits are set:
POPCOUNT8 = [0] * 2**8
for index in xrange(len(POPCOUNT8)):
POPCOUNT8[index] = (index & 1) + POPCOUNT8[index >> 1]
def popcount32(x):
return POPCOUNT8[x & 0xff] + POPCOUNT8[(x >> 8) & 0xff] + POPCOUNT8[(x >> 16) & 0xff] + POPCOUNT8[(x >> 24) & 0xff]
Spin-$1/2$ basis for given number of up-spins.
INDEX_TYPE = np.int64
INVALID_INDEX = -1
class SpinOneHalfBasis:
def __init__(self, num_spins, num_ups):
assert 0 <= num_ups <= num_spins
assert num_spins <= 32
self.num_spins = num_spins
self.num_ups = num_ups
# maps state to its index in the state array
self.states = []
self.index = INVALID_INDEX * np.ones(2**num_spins, dtype=INDEX_TYPE)
for s in xrange(2**num_spins):
if popcount32(s) == num_ups:
self.index[s] = len(self.states)
self.states.append(s)
self.dim = len(self.states)
Heisenberg Hamiltonian:
def heisenberg_hamiltonian(num_spins, num_ups, J, periodic=False):
basis = SpinOneHalfBasis(num_spins, num_ups)
num_pairs = num_spins - 1
D = basis.dim
def matvec(v):
w = np.empty(D)
# \sum_j sigma_z^j sigma_z^(j+1)
mask = 2**num_pairs - 1
for idx, s in enumerate(basis.states):
# count number of parallel & antiparallel pairs of spins
p = popcount32(mask & ~(s ^ (s>>1)))
a = num_pairs - p
# set diagonal entry
w[idx] = 0.25 * (p - a) * v[idx]
# \sum_j sigma_+^j sigma_-^(j+1) + h.c.
for idx, s in enumerate(basis.states):
for r in xrange(num_pairs):
# flip [01, 10] -> 0.5 * [10, 01]
idx2 = basis.index[s ^ (0b11 << r)] # will be INVALID_INDEX if we had [00, 11] at this position, since then the number of up-spins is changed
if idx2 != INVALID_INDEX:
w[idx2] += 0.5 * v[idx]
# sigma_+^1 sigma_-^L + h.c.
if periodic:
for idx, s in enumerate(basis.states):
# diagonal
p = 1 & ~(s ^ (s>>num_pairs))
a = 1 - p
w[idx] += 0.25 * (p - a) * v[idx]
# offdiagonal
idx2 = basis.index[s ^ (1 << num_pairs | 1)]
if idx2 != INVALID_INDEX:
w[idx2] += 0.5 * v[idx]
return J * w
return la.LinearOperator((D, D), matvec, dtype=np.float64)
assert np.isclose(heisenberg_hamiltonian(2, 2, J=1) * [1], [0.25])
assert np.isclose(heisenberg_hamiltonian(2, 0, J=1) * [1], [0.25])
assert np.allclose(heisenberg_hamiltonian(2, 1, J=1) * [1, 0], [-0.25, 0.5])
assert np.allclose(heisenberg_hamiltonian(2, 1, J=1) * [0, 1], [0.5, -0.25])
Compute energies of lowest eigenstates:
def heisenberg_energies(k, num_spins, num_ups, J, periodic):
H = heisenberg_hamiltonian(num_spins, num_ups, J, periodic)
# 1x1 Hamiltonian
if H.shape == (1, 1):
return H * [1] / num_spins
# compute using lanczos
Es = la.eigsh(H, k, which='SA', return_eigenvectors=False)
return Es / num_spins
heisenberg_energies(1, num_spins=2, num_ups=1, J=1, periodic=False)
%%time
for L in range(2, 10):
print L, heisenberg_energies(1, L, L/2, 1, periodic=False)
%%time
for L in range(2, 10):
print L, heisenberg_energies(1, L, L/2, 1, periodic=True)