1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
|
# http://scitation.aip.org/content/aapt/journal/ajp/83/10/10.1119/1.4923249
# Asymmetric wave functions from tiny perturbations
'''
Abstract:
The quantum mechanical behavior of a particle in a double well defies
our intuition based on classical reasoning. Not surprisingly, an asymmetry
in the double well will restore results more consistent with the
classical picture. What is surprising, however, is how a very small asymmetry
can lead to essentially classical behavior. In this paper, we use the simplest
version of a double-well potential to demonstrate these statements.
We also show how this system accurately maps onto a two-state system,
which we refer to as a 'toy model.
License: GPL-V3
Author-email; jithinbp@gmail.com
'''
from scipy import linalg,integrate
import numpy as np
import time
from pylab import *
import pyperclip
def lazycopy(A):
np.savetxt('tmp.csv',A)
fo = open('tmp.csv', 'r').read()
pyperclip.copy(fo)
########## PARAMETERS
spatial_scaling = 1.#0.52917725
a = 1. # width of the infinite 1D potential
b = .2 #width of the barrier
N = 200 #the number of terms in the basis fn expansion
q = 5 #the number of wavefns to be plotted
a = a/spatial_scaling
b = b/spatial_scaling
E10 = np.pi*np.pi/2./a #4.9348022 #(pi^2)/(2*a^2)
print ('E10 = ',E10)
energy_scaling = 1./E10 #27.211396
V0 = 500 #depth of the well
V0= V0/energy_scaling
n = 5000
x = np.linspace(0,a,n)
V = V0*np.zeros(n);
b1= int(n*(a/2 - b/2)/a )
b2= int(n*(a/2 + b/2)/a )
V[b1:b2]=V0;
h=np.zeros([N,N]); #initialising the h matrix
areasR = []
areasL = []
######## PREPARE PLOTS
########Plot ground state + potential's shape
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
ax1.set_xlabel('x/a')
ax1.set_ylabel('psi^2', color='g')
ax2.set_ylabel('V', color='b')
ax2.plot(x*spatial_scaling,V*energy_scaling, 'b-')
ax2.set_ylim([0,V0*energy_scaling+10])
#plt.title(lbl)
############# asymmetry ###########
#delta = np.linspace(0,0.05,6)
delta = [0,2e-7,5e-7,1e-6,2e-6,1e-5]#np.linspace(0,0.05,6)
for VL in delta:
V[0:b1]= -1*VL/energy_scaling
#Determing the elements of the h matrix
for i in range(N):
for j in range(i,N):
bi = np.sqrt(2/a)*np.sin(x*np.float((i+1)*np.pi/a));
bj = np.sqrt(2/a)*np.sin(x*np.float((j+1)*np.pi/a));
f = bi*V*bj;
h[j,i] = integrate.simps(f,x);
if (j==i):
h[j,i] = h[j,i]+((j+1)**2)*(np.pi**2)*(0.5/a**2);
h[i,j]=h[j,i];
evals,evecs = linalg.eig(h)
sortorder = evals.argsort()
E = evals[sortorder]
E=E*energy_scaling;
evecs=evecs[:,sortorder] #Corresonding eigen vectors
#lazycopy(np.diagonal(h,1))
psi = np.zeros([N,len(x)])
for j in range(1): #We only need the ground state for now
for m in range(N):
psi[j,:]=psi[j,:]+evecs[m,j]*np.sqrt(2/a)*np.sin((m+1)*np.pi*x/a) #wavefn using eigen expansion
PY = psi[0,:]**2 # +E[0]
lbl = 'VL = -%.2e'%(VL)
print min(V),E[0],E[1],max(PY[:len(PY)/2]),max(PY[len(PY)/2:])
if VL==0:
ax1.plot(x[::10]*spatial_scaling,PY[::10],color='R',label='symmetric wells')
else:
ax1.plot(x[::10]*spatial_scaling,PY[::10],color='G',label=lbl)
areasL.append(np.trapz(PY[:len(PY)/2],dx = (x[1]-x[0]) ))
areasR.append(np.trapz(PY[len(PY)/2:],dx = (x[1]-x[0]) ))
legend = ax1.legend(loc='upper right', shadow=True)
# The frame is matplotlib.patches.Rectangle instance surrounding the legend.
frame = legend.get_frame()
frame.set_facecolor('0.90')
# Set the fontsize
for label in legend.get_texts():
label.set_fontsize('small')
figure()
xlabel('perturbation(VR-VL)')
ylabel('Area under the curve')
plot(delta,areasL,color='R',label='Left well area')
plot(delta,areasR,color='G',label='Right well area')
show()
|