|
| 1 | +"""Lanczos algortithm with exact diagonalizatoin. |
| 2 | + |
| 3 | +Also: generate Hamiltonian of the transverse field Ising model. |
| 4 | + |
| 5 | +H = -J sum_i sigma^x_i sigma^x_{i+1} - g sum_i sigma^z i; periodic boundary cond. |
| 6 | + |
| 7 | +""" |
| 8 | + |
| 9 | +import numpy as np |
| 10 | +import scipy |
| 11 | +from scipy import sparse |
| 12 | +import scipy.sparse.linalg |
| 13 | +import matplotlib.pyplot as plt |
| 14 | + |
| 15 | +Id = sparse.csr_matrix(np.eye(2)) |
| 16 | +Sx = sparse.csr_matrix([[0., 1.], [1., 0.]]) |
| 17 | +Sz = sparse.csr_matrix([[1., 0.], [0., -1.]]) |
| 18 | +Splus = sparse.csr_matrix([[0., 1.], [0., 0.]]) |
| 19 | +Sminus = sparse.csr_matrix([[0., 0.], [1., 0.]]) |
| 20 | + |
| 21 | + |
| 22 | +def singesite_to_full(op, i, L): |
| 23 | + op_list = [Id]*L # = [Id, Id, Id ...] with L entries |
| 24 | + op_list[i] = op |
| 25 | + full = op_list[0] |
| 26 | + for op_i in op_list[1:]: |
| 27 | + full = sparse.kron(full, op_i, format="csr") |
| 28 | + return full |
| 29 | + |
| 30 | + |
| 31 | +def gen_sx_list(L): |
| 32 | + return [singesite_to_full(Sx, i, L) for i in range(L)] |
| 33 | + |
| 34 | + |
| 35 | +def gen_sz_list(L): |
| 36 | + return [singesite_to_full(Sz, i, L) for i in range(L)] |
| 37 | + |
| 38 | + |
| 39 | +def gen_hamiltonian(sx_list, sz_list, g, J=1.): |
| 40 | + L = len(sx_list) |
| 41 | + H = sparse.csr_matrix((2**L, 2**L)) |
| 42 | + for j in range(L): |
| 43 | + H = H - J *( sx_list[j] * sx_list[(j+1)%L]) |
| 44 | + H = H - g * sz_list[j] |
| 45 | + return H |
| 46 | + |
| 47 | + |
| 48 | +def lanczos(psi0, H, N=200, stabilize=False): |
| 49 | + """Perform a Lanczos iteration building the tridiagonal matrix T and ONB of the Krylov space.""" |
| 50 | + if psi0.ndim != 1: |
| 51 | + raise ValueError("psi0 should be a vector, " |
| 52 | + "i.e., a numpy array with a single dimension of len 2**L") |
| 53 | + if H.shape[1] != psi0.shape[0]: |
| 54 | + raise ValueError("Shape of H doesn't match len of psi0.") |
| 55 | + psi0 = psi0/np.linalg.norm(psi0) |
| 56 | + vecs = [psi0] |
| 57 | + T = np.zeros((N, N)) |
| 58 | + psi = H @ psi0 # @ means matrix multiplication |
| 59 | + # and works both for numpy arrays and scipy.sparse.csr_matrix |
| 60 | + alpha = T[0, 0] = np.inner(psi0.conj(), psi).real |
| 61 | + psi = psi - alpha* vecs[-1] |
| 62 | + for i in range(1, N): |
| 63 | + beta = np.linalg.norm(psi) |
| 64 | + if beta < 1.e-13: |
| 65 | + print("Lanczos terminated early after i={i:d} steps:" |
| 66 | + "full Krylov space built".format(i=i)) |
| 67 | + T = T[:i, :i] |
| 68 | + break |
| 69 | + psi /= beta |
| 70 | + # note: mathematically, psi should be orthogonal to all other states in `vecs` |
| 71 | + if stabilize: |
| 72 | + for vec in vecs: |
| 73 | + psi -= vec * np.inner(vec.conj(), psi) |
| 74 | + psi /= np.linalg.norm(psi) |
| 75 | + vecs.append(psi) |
| 76 | + psi = H @ psi - beta * vecs[-2] |
| 77 | + alpha = np.inner(vecs[-1].conj(), psi).real |
| 78 | + psi = psi - alpha * vecs[-1] |
| 79 | + T[i, i] = alpha |
| 80 | + T[i-1, i] = T[i, i-1] = beta |
| 81 | + return T, vecs |
| 82 | + |
| 83 | + |
| 84 | +def colorplot(xs, ys, data, **kwargs): |
| 85 | + """Create a colorplot with matplotlib.pyplot.imshow. |
| 86 | + |
| 87 | + Parameters |
| 88 | + ---------- |
| 89 | + xs : 1D array, shape (n,) |
| 90 | + x-values of the points for which we have data; evenly spaced |
| 91 | + ys : 1D array, shape (m,) |
| 92 | + y-values of the points for which we have data; evenly spaced |
| 93 | + data : 2D array, shape (m, n) |
| 94 | + ``data[i, j]`` corresponds to the points ``(xs[i], ys[j])`` |
| 95 | + **kwargs : |
| 96 | + additional keyword arguments, given to `imshow`. |
| 97 | + """ |
| 98 | + data = np.asarray(data) |
| 99 | + if data.shape != (len(xs), len(ys)): |
| 100 | + raise ValueError("Shape of data doesn't match len of xs and ys!") |
| 101 | + dx = (xs[-1] - xs[0])/(len(xs)-1) |
| 102 | + assert abs(dx - (xs[1]-xs[0])) < 1.e-10 |
| 103 | + dy = (ys[-1] - ys[0])/(len(ys)-1) |
| 104 | + assert abs(dy - (ys[1]-ys[0])) < 1.e-10 |
| 105 | + extent = (xs[0] - 0.5 * dx, xs[-1] + 0.5 * dx, # left, right |
| 106 | + ys[0] - 0.5 * dy, ys[-1] + 0.5 * dy) # bottom, top |
| 107 | + kwargs.setdefault('aspect', 'auto') |
| 108 | + kwargs.setdefault('interpolation', 'nearest') |
| 109 | + kwargs.setdefault('extent', extent) |
| 110 | + # convention of imshow: matrix like data[row, col] with (0, 0) top left. |
| 111 | + # but we want data[col, row] with (0, 0) bottom left -> transpose and invert y axis |
| 112 | + plt.imshow(data.T[::-1, :], **kwargs) |
0 commit comments