Source code for quimb.gen.operators

"""Functions for generating quantum operators.
"""
from operator import add
import math
import functools
import itertools
import operator

import numpy as np
import scipy.sparse as sp
from scipy.special import comb

from ..utils import isiterable, concat, unique
from ..core import (qarray, make_immutable, get_thread_pool,
                    par_reduce, isreal, qu, eye, kron, ikron)


# --------------------------------------------------------------------------- #
#                      gates and other simple operators                       #
# --------------------------------------------------------------------------- #

[docs]@functools.lru_cache(maxsize=16) def spin_operator(label, S=1 / 2, **kwargs): """Generate a general spin-operator. Parameters ---------- label : str The type of operator, can be one of six options: - ``{'x', 'X'}``, x-spin operator. - ``{'y', 'Y'}``, y-spin operator. - ``{'z', 'Z'}``, z-spin operator. - ``{'+', 'p'}``, Raising operator. - ``{'-', 'm'}``, Lowering operator. - ``{'i', 'I'}``, identity operator. S : float, optional The spin of particle to act on, default to spin-1/2. kwargs Passed to :func:`quimbify`. Returns ------- S : immutable operator The spin operator. See Also -------- pauli Examples -------- >>> spin_operator('x') qarray([[0. +0.j, 0.5+0.j], [0.5+0.j, 0. +0.j]]) >>> qu.spin_operator('+', S=1) qarray([[0. +0.j, 1.41421356+0.j, 0. +0.j], [0. +0.j, 0. +0.j, 1.41421356+0.j], [0. +0.j, 0. +0.j, 0. +0.j]]) >>> qu.spin_operator('Y', sparse=True) <2x2 sparse matrix of type '<class 'numpy.complex128'>' with 2 stored elements in Compressed Sparse Row format> """ D = int(2 * S + 1) op = np.zeros((D, D), dtype=complex) ms = np.linspace(S, -S, D) label = label.lower() if label in {'x', 'y'}: for i in range(D - 1): c = 0.5 * (S * (S + 1) - (ms[i] * ms[i + 1]))**0.5 op[i, i + 1] = -1.0j * c if (label == 'y') else c op[i + 1, i] = 1.0j * c if (label == 'y') else c elif label == 'z': for i in range(D): op[i, i] = ms[i] elif label in {'+', 'p', '-', 'm'}: for i in range(D - 1): c = (S * (S + 1) - (ms[i] * ms[i + 1]))**0.5 if label in {'+', 'p'}: op[i, i + 1] = c else: op[i + 1, i] = c elif label in {'i', 'I'}: np.fill_diagonal(op, 1.0) else: raise ValueError(f"Label '{label}'' not understood, should be one of " "``['X', 'Y', 'Z', '+', '-', 'I']``.") op = qu(np.real_if_close(op), **kwargs) make_immutable(op) return op
[docs]@functools.lru_cache(maxsize=8) def pauli(xyz, dim=2, **kwargs): """Generates the pauli operators for dimension 2 or 3. Parameters ---------- xyz : str Which spatial direction, upper or lower case from ``{'I', 'X', 'Y', 'Z'}``. dim : int, optional Dimension of spin operator (e.g. 3 for spin-1), defaults to 2 for spin half. kwargs Passed to ``quimbify``. Returns ------- P : immutable operator The pauli operator. See Also -------- spin_operator """ xyzmap = {0: 'i', 'i': 'i', 'I': 'i', 1: 'x', 'x': 'x', 'X': 'x', 2: 'y', 'y': 'y', 'Y': 'y', 3: 'z', 'z': 'z', 'Z': 'z'} opmap = {('i', 2): lambda: eye(2, **kwargs), ('x', 2): lambda: qu([[0, 1], [1, 0]], **kwargs), ('y', 2): lambda: qu([[0, -1j], [1j, 0]], **kwargs), ('z', 2): lambda: qu([[1, 0], [0, -1]], **kwargs), ('i', 3): lambda: eye(3, **kwargs), ('x', 3): lambda: qu([[0, 1, 0], [1, 0, 1], [0, 1, 0]], **kwargs) / 2**.5, ('y', 3): lambda: qu([[0, -1j, 0], [1j, 0, -1j], [0, 1j, 0]], **kwargs) / 2**.5, ('z', 3): lambda: qu([[1, 0, 0], [0, 0, 0], [0, 0, -1]], **kwargs)} op = opmap[(xyzmap[xyz], dim)]() # Operator is cached, so make sure it cannot be modified make_immutable(op) return op
[docs]@functools.lru_cache(8) def hadamard(dtype=complex, sparse=False): """The Hadamard gate. """ H = qu([[1., 1.], [1., -1.]], dtype=dtype, sparse=sparse) / 2**0.5 make_immutable(H) return H
[docs]@functools.lru_cache(128) def phase_gate(phi, dtype=complex, sparse=False): """The generalized qubit phase-gate, which adds phase ``phi`` to the ``|1>`` state. """ Rp = qu([[1., 0.], [0., np.exp(1.0j * phi)]], dtype=dtype, sparse=sparse) make_immutable(Rp) return Rp
[docs]@functools.lru_cache(8) def T_gate(dtype=complex, sparse=False): """The T-gate (pi/8 gate). """ return phase_gate(math.pi / 4, dtype=dtype, sparse=sparse)
[docs]@functools.lru_cache(8) def S_gate(dtype=complex, sparse=False): """The S-gate (phase gate). """ return phase_gate(math.pi / 2, dtype=dtype, sparse=sparse)
[docs]@functools.lru_cache(128) def rotation(phi, xyz='Z', dtype=complex, sparse=False): """The single qubit rotation gate. """ R = math.cos(phi / 2) * pauli('I') - 1.0j * math.sin(phi / 2) * pauli(xyz) R = qu(R, dtype=dtype, sparse=sparse) make_immutable(R) return R
Rx = functools.partial(rotation, xyz='x') Ry = functools.partial(rotation, xyz='y') Rz = functools.partial(rotation, xyz='z')
[docs]@functools.lru_cache(128) def U_gate(theta, phi, lamda, dtype=complex, sparse=False): r"""Arbitrary unitary single qubit gate. .. math:: U_3(\theta, \phi, \lambda) = \begin{bmatrix} \cos(\theta / 2) & - e^{i \lambda} \sin(\theta / 2) \\ e^{i \phi} \sin(\theta / 2) & e^{i(\lambda + \phi)}\cos(\theta / 2) \end{bmatrix} Parameters ---------- theta : float Angle between 0 and pi. phi : float Angle between 0 and 2 pi. lamba : float Angle between 0 and 2 pi. Returns ------- U : (2, 2) array The unitary matrix, cached. """ from cmath import cos, sin, exp c2, s2 = cos(theta / 2), sin(theta / 2) U = qu( [[c2, -exp(1j * lamda) * s2], [exp(1j * phi) * s2, exp(1j * (lamda + phi)) * c2]], dtype=dtype, sparse=sparse ) make_immutable(U) return U
[docs]@functools.lru_cache(4) def Xsqrt(**qu_opts): r"""Rx(pi / 2). .. math:: X^{\frac{1}{2}} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & - i \\ - i & 1 \end{bmatrix} """ return Rx(math.pi / 2, **qu_opts)
[docs]@functools.lru_cache(4) def Ysqrt(**qu_opts): r"""Ry(pi / 2). .. math:: Y^{\frac{1}{2}} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & - 1 \\ 1 & 1 \end{bmatrix} """ return Ry(math.pi / 2, **qu_opts)
[docs]@functools.lru_cache(4) def Zsqrt(**qu_opts): r"""Rz(pi / 2). .. math:: Z^{\frac{1}{2}} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 - i & 0 \\ 0 & 1 + i \end{bmatrix} """ return Ry(math.pi / 2, **qu_opts)
[docs]@functools.lru_cache(4) def Wsqrt(**qu_opts): r"""R[X + Y](pi / 2). .. math:: W^{\frac{1}{2}} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & -\sqrt{i} \\ \sqrt{-i} & 1 \end{bmatrix} """ return U_gate(math.pi / 2, -math.pi / 4, math.pi / 4, **qu_opts)
[docs]@functools.lru_cache(maxsize=8) def swap(dim=2, dtype=complex, **kwargs): """The SWAP operator acting on subsystems of dimension `dim`. """ S = np.identity(dim**2, dtype=dtype) S = (S.reshape([dim, dim, dim, dim]) .transpose([0, 3, 1, 2]) .reshape([dim**2, dim**2])) S = qu(S, dtype=dtype, **kwargs) make_immutable(S) return S
[docs]@functools.lru_cache(maxsize=128) def fsim(theta, phi, dtype=complex, **kwargs): r"""The 'fermionic simulation' gate: .. math:: \mathrm{fsim}(\theta, \phi) = \begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & \cos(\theta) & -i sin(\theta) & 0\\ 0 & -i sin(\theta) & \cos(\theta) & 0\\ 0 & 0 & 0 & \exp(-i \phi) \end{bmatrix} Note that ``theta`` and ``phi`` should be specified in radians and the sign convention with this gate varies. Here for example, ``fsim(- pi / 2, 0) == iswap()``. """ from cmath import cos, sin, exp a = cos(theta) b = -1j * sin(theta) c = exp(-1j * phi) gate = [[1, 0, 0, 0], [0, a, b, 0], [0, b, a, 0], [0, 0, 0, c]] gate = qu(gate, dtype=dtype, **kwargs) make_immutable(gate) return gate
@functools.lru_cache(maxsize=4) def iswap(dtype=complex, **kwargs): iswap = qu([[1., 0., 0., 0.], [0., 0., 1j, 0.], [0., 1j, 0., 0.], [0., 0., 0., 1.]], dtype=dtype, **kwargs) make_immutable(iswap) return iswap
[docs]@functools.lru_cache(maxsize=16) def controlled(s, dtype=complex, sparse=False): """Construct a controlled pauli gate for two qubits. Parameters ---------- s : str Which pauli to use, including 'not' aliased to 'x'. sparse : bool, optional Whether to construct a sparse operator. Returns ------- C : immutable operator The controlled two-qubit gate operator. """ # alias not and NOT to x s = {'NOT': 'x', 'not': 'x'}.get(s, s) kws = {'dtype': dtype, 'sparse': sparse} op = ((qu([1, 0], qtype='dop', **kws) & eye(2, **kws)) + (qu([0, 1], qtype='dop', **kws) & pauli(s, **kws))) make_immutable(op) return op
[docs]@functools.lru_cache(8) def CNOT(dtype=complex, sparse=False): """The controlled-not gate. """ return controlled('not', dtype=dtype, sparse=sparse)
[docs]@functools.lru_cache(8) def cX(dtype=complex, sparse=False): """The controlled-X gate. """ return controlled('not', dtype=dtype, sparse=sparse)
[docs]@functools.lru_cache(8) def cY(dtype=complex, sparse=False): """The controlled-Y gate. """ return controlled('Y', dtype=dtype, sparse=sparse)
[docs]@functools.lru_cache(8) def cZ(dtype=complex, sparse=False): """The controlled-Z gate. """ return controlled('Z', dtype=dtype, sparse=sparse)
# --------------------------------------------------------------------------- # # Hamiltonians # # --------------------------------------------------------------------------- #
[docs]def hamiltonian_builder(fn): """Wrap a function to perform some generic postprocessing and take the kwargs ``stype`` and ``sparse``. This assumes the core function always builds the hamiltonian in sparse form. The wrapper then: 1. Checks if the operator is real and, if so, discards imaginary part if no explicity `dtype` was given 2. Converts the operator to dense or the correct sparse form 3. Makes the operator immutable so it can be safely cached """ @functools.wraps(fn) def ham_fn(*args, stype='csr', sparse=False, **kwargs): H = fn(*args, **kwargs) if kwargs.get('dtype', None) is None and isreal(H): H = H.real if not sparse: H = qarray(H.A) elif H.format != stype: H = H.asformat(stype) make_immutable(H) return H return ham_fn
[docs]@functools.lru_cache(maxsize=8) @hamiltonian_builder def ham_heis(n, j=1.0, b=0.0, cyclic=False, parallel=False, nthreads=None, ownership=None): """Constructs the nearest neighbour 1d heisenberg spin-1/2 hamiltonian. Parameters ---------- n : int Number of spins. j : float or tuple(float, float, float), optional Coupling constant(s), with convention that positive = antiferromagnetic. Can supply scalar for isotropic coupling or vector ``(jx, jy, jz)``. b : float or tuple(float, float, float), optional Magnetic field, defaults to z-direction only if tuple not given. cyclic : bool, optional Whether to couple the first and last spins. sparse : bool, optional Whether to return the hamiltonian in sparse form. stype : str, optional What format of sparse operator to return if ``sparse``. parallel : bool, optional Whether to build the operator in parallel. By default will do this for n > 16. nthreads : int optional How mny threads to use in parallel to build the operator. ownership : (int, int), optional If given, which range of rows to generate. kwargs Supplied to :func:`~quimb.core.quimbify`. Returns ------- H : immutable operator The Hamiltonian. """ dims = (2,) * n try: jx, jy, jz = j except TypeError: jx = jy = jz = j try: bx, by, bz = b except TypeError: bz = b bx = by = 0.0 parallel = (n > 16) if parallel is None else parallel op_kws = {'sparse': True, 'stype': 'coo'} ikron_kws = {'sparse': True, 'stype': 'coo', 'coo_build': True, 'ownership': ownership} # The basic operator (interaction and single b-field) that can be repeated. two_site_term = sum( j * kron(spin_operator(s, **op_kws), spin_operator(s, **op_kws)) for j, s in zip((jx, jy, jz), 'xyz') ) - sum( b * kron(spin_operator(s, **op_kws), eye(2, **op_kws)) for b, s in zip((bx, by, bz), 'xyz') if b != 0.0 ) single_site_b = sum(-b * spin_operator(s, **op_kws) for b, s in zip((bx, by, bz), 'xyz') if b != 0.0) def gen_term(i): # special case: the last b term needs to be added manually if i == -1: return ikron(single_site_b, dims, n - 1, **ikron_kws) # special case: the interaction between first and last spins if cyclic if i == n - 1: return sum( j * ikron(spin_operator(s, **op_kws), dims, [0, n - 1], **ikron_kws) for j, s in zip((jx, jy, jz), 'xyz') if j != 0.0) # General term, on-site b-field plus interaction with next site return ikron(two_site_term, dims, [i, i + 1], **ikron_kws) terms_needed = range(0 if not any((bx, by, bz)) else -1, n if cyclic else n - 1) if parallel: pool = get_thread_pool(nthreads) ham = par_reduce(add, pool.map(gen_term, terms_needed)) else: ham = sum(map(gen_term, terms_needed)) return ham
[docs]def ham_ising(n, jz=1.0, bx=1.0, **ham_opts): """Generate the quantum transverse field ising model hamiltonian. This is a simple alias for :func:`~quimb.gen.operators.ham_heis` with Z-interactions and an X-field. """ return ham_heis(n, j=(0, 0, jz), b=(bx, 0, 0), **ham_opts)
[docs]def ham_XY(n, jxy, bz, **ham_opts): """Generate the quantum transverse field XY model hamiltonian. This is a simple alias for :func:`~quimb.gen.operators.ham_heis` with X- and Y-interactions and a Z-field. """ return ham_heis(n, j=(jxy, jxy, 0), b=(0, 0, bz), **ham_opts)
[docs]def ham_XXZ(n, delta, jxy=1.0, **ham_opts): """Generate the XXZ-model hamiltonian. This is a simple alias for :func:`~quimb.gen.operators.ham_heis` with matched X- and Y-interactions and ``delta`` Z coupling. """ return ham_heis(n, j=(jxy, jxy, delta), b=0, **ham_opts)
[docs]@functools.lru_cache(maxsize=8) @hamiltonian_builder def ham_j1j2(n, j1=1.0, j2=0.5, bz=0.0, cyclic=False, ownership=None): """Generate the j1-j2 hamiltonian, i.e. next nearest neighbour interactions. Parameters ---------- n : int Number of spins. j1 : float, optional Nearest neighbour coupling strength. j2 : float, optional Next nearest neighbour coupling strength. bz : float, optional B-field strength in z-direction. cyclic : bool, optional Cyclic boundary conditions. sparse : bool, optional Return hamiltonian as sparse-csr operator. ownership : (int, int), optional If given, which range of rows to generate. kwargs Supplied to :func:`~quimb.core.quimbify`. Returns ------- H : immutable operator The Hamiltonian. """ dims = (2,) * n op_kws = {'sparse': True, 'stype': 'coo'} ikron_kws = {'sparse': True, 'stype': 'coo', 'coo_build': True, 'ownership': ownership} sxyz = [spin_operator(i, **op_kws) for i in 'xyz'] coosj1 = np.array([(i, i + 1) for i in range(n)]) coosj2 = np.array([(i, i + 2) for i in range(n)]) if cyclic: coosj1, coosj2 = coosj1 % n, coosj2 % n else: coosj1 = coosj1[np.all(coosj1 < n, axis=1)] coosj2 = coosj2[np.all(coosj2 < n, axis=1)] def j1_terms(): for coo in coosj1: if abs(coo[1] - coo[0]) == 1: # can sum then tensor (faster) yield ikron(sum(op & op for op in sxyz), dims, coo, **ikron_kws) else: # tensor then sum (slower) yield sum(ikron(op, dims, coo, **ikron_kws) for op in sxyz) def j2_terms(): for coo in coosj2: if abs(coo[1] - coo[0]) == 2: # can add then tensor (faster) yield ikron(sum(op & eye(2, **op_kws) & op for op in sxyz), dims, coo, **ikron_kws) else: yield sum(ikron(op, dims, coo, **ikron_kws) for op in sxyz) ham = j1 * sum(j1_terms()) + j2 * sum(j2_terms()) if bz != 0: gen_bz = (ikron([sxyz[2]], dims, i, **ikron_kws) for i in range(n)) ham += bz * sum(gen_bz) return ham
def _gen_mbl_random_factors(n, dh, dh_dim, dh_dist, seed=None, beta=None): # sort out a vector of noise strengths -> e.g. (0, 0, 1) for z-noise only if isinstance(dh, (tuple, list)): dhds = dh else: dh_dim = {0: '', 1: 'z', 2: 'xy', 3: 'xyz'}.get(dh_dim, dh_dim) dhds = tuple((dh if d in dh_dim else 0) for d in 'xyz') if seed is not None: np.random.seed(seed) # sort out the noise distribution if dh_dist in {'g', 'gauss', 'gaussian', 'normal'}: rs = np.random.randn(3, n) elif dh_dist in {'s', 'flat', 'square', 'uniform', 'box'}: rs = 2.0 * np.random.rand(3, n) - 1.0 elif dh_dist in {'qp', 'quasiperiodic', 'qr', 'quasirandom'}: if dh_dim != 'z': raise ValueError("dh_dim should be 1 or 'z' for dh_dist='qp'.") if beta is None: beta = (5**0.5 - 1) / 2 # the random phase delta = 2 * np.pi * np.random.rand() # make sure get 3 by n different strengths inds = np.broadcast_to(range(n), (3, n)) rs = np.cos(2 * np.pi * beta * inds + delta) return dhds, rs
[docs]@hamiltonian_builder def ham_mbl(n, dh, j=1.0, bz=0.0, cyclic=False, seed=None, dh_dist="s", dh_dim=1, beta=None, ownership=None): """ Constructs a heisenberg hamiltonian with isotropic coupling and random fields acting on each spin - the many-body localized (MBL) spin hamiltonian. Parameters ---------- n : int Number of spins. dh : float or (float, float, float) Strength of random fields (stdev of gaussian distribution), can be scalar (isotropic noise) or 3-vector for (x, y, z) directions. j : float or (float, float, float), optional Coupling strength, can be scalar (isotropic) or 3-vector. bz : float, optional Global magnetic field (in z-direction). cyclic : bool, optional Whether to use periodic boundary conditions. seed : int, optional Number to seed random number generator with. dh_dist : {'g', 's', 'qr'}, optional Type of random distribution for the noise: - "s": square, with bounds ``(-dh, dh)`` - "g": gaussian, with standard deviation ``dh`` - "qp": quasi periodic, with amplitude ``dh`` and 'wavenumber' ``beta`` so that the field at site ``i`` is ``dh * cos(2 * pi * beta * i + delta)`` with ``delta`` a random offset between ``(0, 2 * pi)``, possibly seeded by ``seed``. dh_dim : {1, 2, 3} or str, optional The number of dimensions the noise acts in, or string specifier like ``'yz'``. beta : float, optional The wave number if ``dh_dist='qr'``, defaults to the golden ratio``(5**0.5 - 1) / 2``. sparse : bool, optional Whether to construct the hamiltonian in sparse form. stype : {'csr', 'csc', 'coo'}, optional The sparse format. ownership : (int, int), optional If given, which range of rows to generate. kwargs Supplied to :func:`~quimb.core.quimbify`. Returns ------- H : operator The MBL hamiltonian for spin-1/2. See Also -------- MPO_ham_mbl """ dhds, rs = _gen_mbl_random_factors(n, dh, dh_dim, dh_dist, seed, beta) # the base hamiltonian ('csr' is most efficient format to add with) ham = ham_heis(n=n, j=j, b=bz, cyclic=cyclic, sparse=True, stype='csr', ownership=ownership) op_kws = {'sparse': True, 'stype': 'coo'} ikron_kws = {'sparse': True, 'stype': 'coo', 'coo_build': True, 'ownership': ownership} def dh_terms(): for i in range(n): # dhd - the total strength in direction x, y, or z # r - the random strength in direction x, y, or z for site i hdh = sum(dhd * r * spin_operator(s, **op_kws) for dhd, r, s in zip(dhds, rs[:, i], 'xyz')) yield ikron(hdh, (2,) * n, i, **ikron_kws) ham = ham + sum(dh_terms()) return ham
[docs]@hamiltonian_builder def ham_heis_2D(n, m, j=1.0, bz=0.0, cyclic=False, parallel=False, ownership=None): """Construct the 2D spin-1/2 heisenberg model hamiltonian. Parameters ---------- n : int The number of rows. m : int The number of columns. j : float or (float, float, float), optional The coupling strength(s). Isotropic if scalar else if vector ``(Jx, Jy, Jz) = j``. bz : float, optional The z direction magnetic field. cyclic : bool, optional Whether to use periodic boundary conditions. sparse : bool, optional Whether to construct the hamiltonian in sparse form. stype : {'csr', 'csc', 'coo'}, optional The sparse format. parallel : bool, optional Construct the hamiltonian in parallel. Faster but might use more memory. ownership : (int, int), optional If given, which range of rows to generate. kwargs Supplied to :func:`~quimb.core.quimbify`. Returns ------- H : operator The hamiltonian. """ # parse interaction strengths try: jx, jy, jz = j except (TypeError, ValueError): jx = jy = jz = j dims = [[2] * m] * n # shape (n, m) sites = tuple(itertools.product(range(n), range(m))) # generate neighbouring pair coordinates def gen_pairs(): for i, j in sites: above, right = (i + 1) % n, (j + 1) % m # ignore wraparound coordinates if not cyclic if cyclic or above != 0: yield ((i, j), (above, j)) if cyclic or right != 0: yield ((i, j), (i, right)) # build the hamiltonian in sparse 'coo' format always for efficiency op_kws = {'sparse': True, 'stype': 'coo'} ikron_kws = {'sparse': True, 'stype': 'coo', 'coo_build': True, 'ownership': ownership} # generate all pairs of coordinates and directions pairs_ss = tuple(itertools.product(gen_pairs(), 'xyz')) # generate XX, YY and ZZ interaction from # e.g. arg ([(3, 4), (3, 5)], 'z') def interactions(pair_s): pair, s = pair_s Sxyz = spin_operator(s, **op_kws) J = {'x': jx, 'y': jy, 'z': jz}[s] return ikron(J * Sxyz, dims, inds=pair, **ikron_kws) # generate Z field def fields(site): Sz = spin_operator('z', **op_kws) return ikron(bz * Sz, dims, inds=[site], **ikron_kws) if not parallel: # combine all terms all_terms = itertools.chain( map(interactions, pairs_ss), map(fields, sites) if bz != 0.0 else ()) H = sum(all_terms) else: pool = get_thread_pool() all_terms = itertools.chain( pool.map(interactions, pairs_ss), pool.map(fields, sites) if bz != 0.0 else ()) H = par_reduce(add, all_terms) return H
[docs]def uniq_perms(xs): """Generate all the unique permutations of sequence ``xs``. Examples -------- >>> list(uniq_perms('0011')) [('0', '0', '1', '1'), ('0', '1', '0', '1'), ('0', '1', '1', '0'), ('1', '0', '0', '1'), ('1', '0', '1', '0'), ('1', '1', '0', '0')] """ if len(xs) == 1: yield (xs[0],) else: uniq_xs = unique(xs) for first_x in uniq_xs: rem_xs = list(xs) rem_xs.remove(first_x) for sub_perm in uniq_perms(rem_xs): yield (first_x,) + sub_perm
[docs]@functools.lru_cache(maxsize=8) def zspin_projector(n, sz=0, stype="csr", dtype=float): """Construct the projector onto spin-z subpspaces. Parameters ---------- n : int Total size of spin system. sz : float or sequence of floats Spin-z value(s) subspace(s) to find projector for. stype : str Sparse format of the output operator. dtype : {float, complex}, optional The data type of the operator to generate. Returns ------- prj : immutable sparse operator, shape (2**n, D) The (non-square) projector onto the specified subspace(s). The subspace size ``D`` is given by ``n choose (n / 2 + s)`` for each ``s`` specified in ``sz``. Examples -------- >>> zspin_projector(n=2, sz=0).A array([[0., 0.], [1., 0.], [0., 1.], [0., 0.]] Project a 9-spin Heisenberg-Hamiltonian into its spin-1/2 subspace: >>> H = ham_heis(9, sparse=True) >>> H.shape (512, 512) >>> P = zspin_projector(n=9, sz=1 / 2) >>> H0 = P.T @ H @ P >>> H0.shape (126, 126) """ if not isiterable(sz): sz = (sz,) p = 0 all_perms = [] for s in sz: # Number of 'up' spins k = n / 2 + s if not k.is_integer(): raise ValueError(f"{s} is not a valid spin half subspace for {n} " "spins.") k = int(round(k)) # Size of subspace p += comb(n, k, exact=True) # Find all computational basis states with correct number of 0s and 1s base_perm = '0' * (n - k) + '1' * k all_perms += [uniq_perms(base_perm)] # Coordinates cis = tuple(range(p)) # arbitrary basis cjs = tuple(int("".join(perm), 2) for perm in concat(all_perms)) # Construct matrix which projects only on to these basis states prj = sp.coo_matrix((np.ones(p, dtype=dtype), (cjs, cis)), shape=(2**n, p), dtype=dtype) prj = qu(prj, stype=stype, dtype=dtype) make_immutable(prj) return prj
[docs]@functools.lru_cache(8) def create(n=2, **qu_opts): """The creation operator acting on an n-level system. """ data = np.zeros((n, n)) for i in range(n): data[i, i - 1] = i ** 0.5 ap = qu(data, **qu_opts) make_immutable(ap) return ap
[docs]@functools.lru_cache(8) def destroy(n=2, **qu_opts): """The annihilation operator acting on an n-level system. """ am = create(n, **qu_opts).T.copy() make_immutable(am) return am
[docs]@functools.lru_cache(8) def num(n, **qu_opts): """The number operator acting on an n-level system. """ ap, am = create(n, **qu_opts), destroy(n, **qu_opts) an = qu(ap @ am, **qu_opts) make_immutable(an) return an
[docs]@functools.lru_cache(maxsize=8) @hamiltonian_builder def ham_hubbard_hardcore(n, t=0.5, V=1., mu=1., cyclic=False, parallel=False, ownership=None): """Generate the spinless fermion hopping hamiltonian. Parameters ---------- n : int The number of sites. t : float, optional The hopping energy. V : float, optional The interaction energy. mu : float, optional The chemical potential - defaults to half-filling. cyclic : bool, optional Whether to use periodic boundary conditions. parallel : bool, optional Construct the hamiltonian in parallel. Faster but might use more memory. ownership : (int, int), optional If given, which range of rows to generate. kwargs Supplied to :func:`~quimb.core.quimbify`. Returns ------- H : operator The hamiltonian. """ op_kws = {'sparse': True, 'stype': 'coo'} ikron_kws = {'sparse': True, 'stype': 'csr', 'coo_build': True, 'ownership': ownership} cdag, c, cnum = (f(2, **op_kws) for f in (create, destroy, num)) neighbor_term = t * ((cdag & c) + (c & cdag)) + V * (cnum & cnum) dims = [2] * n def terms(): # interacting terms for i, j in [(i, i + 1) for i in range(n - 1)]: yield ikron(neighbor_term, dims, (i, j), **ikron_kws) if cyclic: # can't sum terms and kron later since identity in middle yield ikron([t * cdag, c], dims, (0, n - 1), **ikron_kws) yield ikron([t * c, cdag], dims, (0, n - 1), **ikron_kws) yield ikron([V * cnum, cnum], dims, (0, n - 1), **ikron_kws) # single site terms for i in range(n): yield ikron(-mu * cnum, dims, i, **ikron_kws) if parallel is None: parallel = (n >= 14) if parallel: return par_reduce(operator.add, terms()) return functools.reduce(operator.add, terms())