Robert Johansson (robert@riken.jp)
%pylab inline
Populating the interactive namespace from numpy and matplotlib
from scipy import *
from scipy import optimize
from wavefunction import *
from wavefunction.wavefunction1d import *
h = 6.626e-34
h_ = h/(2*pi)
e = 1.602e-19
Phi0 = h / (2 * e)
cf = 1
ceta = 0.05
J = 0.9725
Ic = 13.3e-6
Ib = J * Ic
Cj = 4.3e-12
mm = Cj * (1 + ceta) * (Phi0/(2*pi))**2 * cf
Ej = Phi0/(2*pi) * Ic / cf
args = {'Ej': Ej, 'Ic': Ic, 'Ib': Ib}
k = -h_ ** 2 / (2 * mm)
x_min = -pi
x_max = pi
N = 750
def U_current_biased(x, args):
"""
Potential for a current-biased phase qubit potential
(the washboard potential)
"""
Ej = args['Ej']
Ic = args['Ic']
Ib = args['Ib']
u = - Ej * (cos(x) + Ib / Ic * x)
return u
x = linspace(x_min, x_max, N)
U = U_current_biased(x, args)
# analytical formula for minima points of potential
x_opt_min = arcsin(Ib/Ic)
print("x_opt_min = %f " % x_opt_min)
x_opt_max = -math.asin(Ib/Ic) + pi
print("x_opt_max = %f" % x_opt_max)
U_min = U_current_biased(x_opt_min, args)
U_max = U_current_biased(x_opt_max, args)
dU = U_max - U_min
dx = x_opt_max - x_opt_min
print("dU = %f" % dU)
x_opt_min = 1.335735 x_opt_max = 1.805858 dU = 0.000000
fig, ax = subplots()
ax.plot(x, U/Ej)
ax.plot(x_opt_min, U_current_biased(x_opt_min, args) / Ej, 'o')
ax.set_ylim(-1.55, -1.5)
ax.set_xlim(0.8,2.2)
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$U(x)/E_J$', fontsize=18);
u = assemble_u_potential(N, U_current_biased, x, args)
K = assemble_K(N, k, x_min, x_max)
V = assemble_V(N, u, x_min, x_max)
H = K + V
evals, evecs = solve_eigenproblem(H)
nbound = 0
boundidx = []
for i in range(0,N):
if evals[i] > U_min - 0.5 * dU and evals[i] < U_max + 0.5 * dU:
if inner(evecs[i] * (x_opt_min-dx < x) * (x < x_opt_max), evecs[i]) > 0.95:
nbound = nbound + 1
boundidx.append(i)
print "Found bound states: ", boundidx
Found bound states: [91, 92, 94, 96, 97, 99, 101, 102, 104]
NN = len(boundidx)
fig, axes = subplots(NN, 1, figsize=(10, NN * 1), sharex=True, sharey=True)
for n, m in enumerate(boundidx):
axes[n].plot(x, real(evecs[m]))
axes[n].set_xlim(0.5, 2.5)
axes[n].set_xlabel(r'$x$', fontsize=18)
fig.tight_layout();
fig, ax = subplots(figsize=(12,8))
ax.plot(x, U/Ej, 'k')
ax.plot(x_opt_min, U_current_biased(x_opt_min, args) / Ej, '.')
for n, m in enumerate(boundidx):
Y0 = evals[m]/Ej * ones(shape(x))
Y = 0.001 * evecs[m] + Y0
ax.plot(x, Y0.real, 'k--')
ax.plot(x, Y.real)
ax.set_ylim(-1.535, -1.52)
ax.set_xlim(0.8,2.2)
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$U(x)$', fontsize=18);
inner_prod = zeros((len(boundidx), len(boundidx))).astype(np.complex)
expect_pos = zeros((len(boundidx), len(boundidx))).astype(np.complex)
expect_kin = zeros((len(boundidx), len(boundidx))).astype(np.complex)
for i, l in enumerate(boundidx):
for j, k in enumerate(boundidx):
psi_l = wavefunction_normalize(x, evecs[l])
psi_k = wavefunction_normalize(x, evecs[k])
inner_prod[i,j] = inner_product(x, psi_l, psi_k) # <psi_l|psi_k>
expect_pos[i,j] = inner_product(x, psi_l, x * psi_k) # <psi_l|x|psi_k>
expect_kin[i,j] = inner_product(x, psi_l, derivative(x, psi_k)) # <psi_l|p|psi_k>
print("Bound energy-levels:")
(evals[boundidx] - evals[boundidx[0]]).real / Ej
Bound energy-levels:
array([ 0. , 0.00107872, 0.00213317, 0.003161 , 0.00415908, 0.005123 , 0.00604603, 0.00691584, 0.00770281])
print("Harmonicity:")
((evals[boundidx] - evals[boundidx[0]]) / (evals[boundidx[1]]-evals[boundidx[0]])).real
Harmonicity:
array([ 0. , 1. , 1.97749688, 2.9303132 , 3.85555056, 4.74913082, 5.60479901, 6.41112924, 7.14067069])
print "inner_prod:\n"
print_matrix(inner_prod)
inner_prod: 1.000 -0.000 -0.000 -0.000 -0.000 0.000 0.000 0.000 -0.000 -0.000 1.000 0.000 -0.000 -0.000 -0.000 -0.000 0.000 0.000 -0.000 0.000 1.000 0.000 -0.000 -0.000 0.000 0.000 0.000 -0.000 -0.000 0.000 1.000 0.000 0.000 -0.000 -0.000 0.000 -0.000 -0.000 -0.000 0.000 1.000 -0.000 0.000 0.000 -0.000 0.000 -0.000 -0.000 0.000 -0.000 1.000 -0.000 0.000 0.000 0.000 -0.000 0.000 -0.000 0.000 -0.000 1.000 -0.000 -0.000 0.000 0.000 0.000 -0.000 0.000 0.000 -0.000 1.000 -0.000 -0.000 0.000 0.000 0.000 -0.000 0.000 -0.000 -0.000 1.000
print "expect_pos:\n"
print_matrix(expect_pos)
expect_pos: 1.341 0.049 -0.002 0.000 -0.000 0.000 -0.000 0.000 -0.000 0.049 1.351 0.070 -0.005 0.000 -0.000 0.000 -0.000 0.000 -0.002 0.070 1.363 0.086 -0.007 0.001 -0.000 0.000 -0.000 0.000 -0.005 0.086 1.376 0.100 -0.010 0.001 -0.000 0.000 -0.000 0.000 -0.007 0.100 1.390 0.113 -0.013 0.002 -0.000 0.000 -0.000 0.001 -0.010 0.113 1.407 0.125 -0.017 0.002 -0.000 0.000 -0.000 0.001 -0.013 0.125 1.428 0.137 -0.023 0.000 -0.000 0.000 -0.000 0.002 -0.017 0.137 1.456 0.146 -0.000 0.000 -0.000 0.000 -0.000 0.002 -0.023 0.146 1.515
print "expect_kin:\n"
print_matrix(expect_kin)
expect_kin: -0.000 10.132 -1.024 0.103 -0.019 0.004 -0.001 0.000 -0.000 -10.132 -0.000 14.099 -1.830 0.215 -0.044 0.012 -0.003 0.001 1.024 -14.099 0.000 16.956 -2.683 0.361 -0.081 0.024 -0.007 -0.103 1.830 -16.956 0.000 19.171 -3.617 0.551 -0.133 0.044 0.019 -0.215 2.683 -19.171 0.000 20.899 -4.670 0.810 -0.207 -0.004 0.044 -0.361 3.617 -20.899 0.000 22.167 -5.909 1.194 0.001 -0.012 0.081 -0.551 4.670 -22.167 -0.000 22.850 -7.450 -0.000 0.003 -0.024 0.133 -0.810 5.909 -22.850 -0.000 22.068 0.000 -0.001 0.007 -0.044 0.207 -1.194 7.450 -22.068 0.000
%reload_ext version_information
%version_information numpy, scipy, matplotlib, wavefunction
Software | Version |
---|---|
Python | 2.7.4 (default, Apr 19 2013, 18:28:01) [GCC 4.7.3] |
IPython | 2.0.0-dev |
OS | posix [linux2] |
numpy | 1.8.0.dev-895866d |
scipy | 0.13.0.dev-7c6d92e |
matplotlib | 1.3.x |
wavefunction | 1.0.0 |
Mon Sep 02 15:16:18 2013 JST |