%matplotlib inline import numpy as np import matplotlib.pyplot as plt import matplotlib.cm as cm from wavefunction.wavefunction2d import * # pick a truncation of the fourier series: n1 = [-L1, ..., L1], n2 = [-L2, ..., L2] L1 = 3 L2 = 3 # pick periods for the coordinates Tx1 = Tx2 = 2 * np.pi # k11, k12, k22 = 0.0000125, 0.0001, 0.0000125 help(assemble_K) K = assemble_K(L1, L2, k11, k12, k22, Tx1, Tx2) args = {'U0': 0, 'eps1': 1, 'eps2': 1, 'eps12': 0.75, 'f': 0.5} x1 = np.linspace(-2*np.pi, 2*np.pi, 100) x2 = np.linspace(-2*np.pi, 2*np.pi, 100) X1, X2 = np.meshgrid(x1, x2) def U_flux_qubit(x1, x2, args): U0, eps1, eps2, eps12, f = args['U0'], args['eps1'], args['eps2'], args['eps12'], args['f'] return U0 - eps1 * np.cos(x1) - eps2 * np.cos(x2) - eps12 * np.cos(2 * np.pi * f + x1 - x2) Z = np.real(U_flux_qubit(X1, X2, args).T) fig, ax = plt.subplots() p = ax.pcolor(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max()) cb = fig.colorbar(p, ax=ax) def assemble_u_flux_qubit(L1, L2, args, sparse=False): U0, eps1, eps2, eps12, f = args['U0'], args['eps1'], args['eps2'], args['eps12'], args['f'] L1n = 2 * L1 + 1 L2n = 2 * L2 + 1 u = np.zeros((L1n, L2n), dtype=np.complex) for n1 in np.arange(-L1, L1+1): N1 = n1 + L1 for n2 in np.arange(-L2, L2+1): N2 = n2 + L2 u[N1, N2] = U0 * delta(0, n1) * delta(0, n2) + \ - 0.5 * eps1 * (delta(1, n1) + delta(-1, n1)) * delta(0, n2) + \ - 0.5 * eps2 * (delta(1, n2) + delta(-1, n2)) * delta(0, n1)+ \ - 0.5 * eps12 * np.exp(+2j * np.pi * f) * delta(+1, n1) * delta(-1, n2) + \ - 0.5 * eps12 * np.exp(-2j * np.pi * f) * delta(-1, n1) * delta(+1, n2) return u def assemble_u_flux_qubit_direct(L1, L2, args, sparse=False): U0, eps1, eps2, eps12, f = args['U0'], args['eps1'], args['eps2'], args['eps12'], args['f'] L1n = 2 * L1 + 1 L2n = 2 * L2 + 1 u = np.zeros((L1n, L2n), dtype=np.complex) u[L1, L2] = U0 u[+1+L1, L2] = - 0.5 * eps1 u[-1+L1, L2] = - 0.5 * eps1 u[L1, +1+L2] = - 0.5 * eps2 u[L1, -1+L2] = - 0.5 * eps2 u[+1+L1, -1+L2] = -0.5 * eps12 * np.exp(+2j * np.pi * f) u[-1+L1, +1+L2] = -0.5 * eps12 * np.exp(-2j * np.pi * f) return np.real(u) u2 = assemble_u_flux_qubit(L1, L2, args) u = assemble_u_flux_qubit_direct(L1, L2, args) abs(u - u2).max() # show code of evalute_fourier_series help(evalute_fourier_series) Z = np.real(evalute_fourier_series(X1, X2, L1, L2, u)) fig, ax = plt.subplots() p = ax.pcolor(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max()) cb = fig.colorbar(p, ax=ax) help(assemble_V) def assemble_V_flux_qubit(L1, L2, args, sparse=False): U0, eps1, eps2, eps12, f = args['U0'], args['eps1'], args['eps2'], args['eps12'], args['f'] L1n = 2 * L1 + 1 L2n = 2 * L2 + 1 V = np.zeros((L1n*L1n, L2n*L2n), dtype=np.complex) for n1 in np.arange(-L1, L1+1): for n2 in np.arange(-L1, L1+1): N = index_m2v(L1, n1, n2) for m1 in np.arange(-L2, L2+1): for m2 in np.arange(-L2, L2+1): M = index_m2v(L2, m1, m2) V[N,M] = U0 * delta(m1, n1) * delta(m2, n2) + \ - 0.5 * eps1 * (delta(m1 + 1, n1) + delta(m1 - 1, n1)) * delta(m2, n2) + \ - 0.5 * eps2 * (delta(m2 + 1, n2) + delta(m2 - 1, n2)) * delta(m1, n1) + \ - 0.5 * eps12 * np.exp(+2j * np.pi * f) * delta(m1 + 1, n1) * delta(m2 - 1, n2) + \ - 0.5 * eps12 * np.exp(-2j * np.pi * f) * delta(m1 - 1, n1) * delta(m2 + 1, n2) return V V_full = assemble_V(L1, L2, u) V = assemble_V_flux_qubit(L1, L2, args) abs(V-V_full).max() H = K + V help(solve_eigenproblem) vals, vecs = solve_eigenproblem(H) fig, ax = plt.subplots() U_x = np.real(U_flux_qubit(x1, -x1, args)) for val in vals: ax.plot(x1, val.real * np.ones_like(x1), 'k') ax.plot(x1, U_x, label="potential", lw='2') ax.axis('tight') ax.set_ylim(min(vals[0], U_x.min()), U_x.max()) ax.set_ylabel(r'$U(x_1, x_2=x_1)$', fontsize=16) ax.set_xlabel(r'$x_1$', fontsize=16); Nstates = 4 fig, axes = plt.subplots(Nstates, 3, figsize=(16, 5 * Nstates), squeeze=False) for n in range(Nstates): psi = convert_v2m(L1, L2, vecs[n]) Z = np.real(evalute_fourier_series(X1, X2, L1, L2, u)) PSI = evalute_fourier_series(X1, X2, L1, L2, psi) p = axes[n, 0].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.real(PSI), cmap=cm.RdBu, vmin=-abs(PSI.real).max(), vmax=abs(PSI.real).max()) c = axes[n, 0].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max()) cb = fig.colorbar(p, ax=axes[n, 0]) cb.set_clim(-5, 5) p = axes[n, 1].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.imag(PSI), cmap=cm.RdBu, vmin=-abs(PSI.imag).max(), vmax=abs(PSI.imag).max()) c = axes[n, 1].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max()) cb = fig.colorbar(p, ax=axes[n, 1]) cb.set_clim(-5, 5) p = axes[n, 2].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.abs(PSI), cmap=cm.Blues, vmin=0, vmax=abs(PSI).max()) c = axes[n, 2].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max()) cb = fig.colorbar(p, ax=axes[n, 2]) axes[1, 0].set_ylabel("%d state" % n); axes[0, 0].set_title(r"$\mathrm{re}\;\Psi(\varphi_p, \varphi_m)$", fontsize=16); axes[0, 1].set_title(r"$\mathrm{im}\;\Psi(\varphi_p, \varphi_m)$", fontsize=16); axes[0, 2].set_title(r"$|\Psi(\varphi_p, \varphi_m)|^2$", fontsize=16); f_vec = np.linspace(0.4, 0.6, 50) e_vals = np.zeros((len(vals), len(f_vec))) for f_idx, f in enumerate(f_vec): K = assemble_K(L1, L2, k11, k12, k22, Tx1, Tx2) args['f'] = f V = assemble_V_flux_qubit(L1, L2, args) H = K + V vals, vecs = solve_eigenproblem(H) e_vals[:, f_idx] = np.real(vals) fig, ax = plt.subplots(1, 1, figsize=(12,6)) for n in range(len(vals)): ax.plot(f_vec, e_vals[n, :]) ax.axis('tight') ax.set_ylim(e_vals[0, :].min(), e_vals[10, :].max()) ax.set_ylabel("Energy levels") ax.set_xlabel(r'$f$'); np.diff(X1, axis=1) np.diff(X2, axis=0).shape psi.shape psi0 = wavefunction_normalize(X1, X2, evalute_fourier_series(X1, X2, L1, L2, convert_v2m(L1, L2, vecs[0]))) psi0.shape derivative(X1, X2, psi0, axis=1) inner_product(X1, X2, psi0, derivative(X1, X2, psi0, axis=0)) Z = np.real(evalute_fourier_series(X1, X2, L1, L2, u)) fig, ax = plt.subplots() p = ax.pcolor(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max()) cb = fig.colorbar(p, ax=ax) Z = X1 ** 2 + X2 ** 2 #Z = np.real(derivative(X1, X2, Z, axis=1)) fig, ax = plt.subplots() p = ax.pcolor(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=-Z.max(), vmax=Z.max()) cb = fig.colorbar(p, ax=ax)