dolfin team mailing list archive

[Chaffra <question142055@xxxxxxxxxxxxxxxxxxxxx>]

New question #142055 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/142055

Dear all,

I am trying to solve the following schrodinger equation in 1D for a quantum well. My problem is that I get different eigenvectors for potential energies that are physically the same. If I use

class BandOffset(Expression):
    def eval(self,v,x):
        v[0] = 0
        if 0.5-0.1  < x[0] < 0.5+0.1: v[0] = -0.23*qe
the wavefunction of eigenstate 0 sits in the well. But if I use

class BandOffset(Expression):
    def eval(self,v,x):
        v[0] = 0.23*qe
        if 0.5-0.1  < x[0] < 0.5+0.1: v[0] = 0
the wavefunction of eigenstate 0 sits in the barrier, but I think the 2 potentials are the same "physically". Thanks for any suggestions. The script is below.



from dolfin import *

kb = 1.3806503e-23
hbar = 1.05457148e-34
qe = 1.60217646e-19
me = 9.10938188e-31

class BandOffset(Expression):
    def eval(self,v,x):
        v[0] = 0
        if 0.5-0.1  < x[0] < 0.5+0.1: v[0] = -0.23*qe

class ElectronEffectiveMass(Expression):
    def eval(self,v,x):
        v[0] = 0.092*me
        if 0.5-0.1 < x[0] < 0.5+0.1: v[0] = 0.067*me
        
def plot_wave(eigensolver, i, Q, **kwargs):
        r, c, rx, cx = eigensolver.get_eigenpair(i)
        print "Eigenvalue ", i,":", r, c

        # Initialize function with eigenvector
        ru = Function(Q, rx)
        cu = Function(Q, cx)
        usq = inner(ru,ru) + inner(cu,cu)
        plot(usq,**kwargs)
        
mesh = Interval(100,0,1)
Q = FunctionSpace(mesh, "CG", 1)

u = TrialFunction(Q)
v = TestFunction(Q)
        
_BO = BandOffset()
_ME = ElectronEffectiveMass()

PE = interpolate(_BO, Q) #potential energy function
ME  = interpolate(_ME, Q) #effective mass function

a = -(hbar**2/2.0)*inner((1.0/ME)*grad(u), grad(v))*dx + PE*v*u*dx
A = PETScMatrix()   
#print "Assembling matrix"
assemble(a, tensor=A) 
#print "Loading eigensolver."
electron = SLEPcEigenSolver()
# Compute all eigenvalues of A x = \lambda x
#print "Computing eigenvalues. This can take a minute."
electron.solve(A)

last_e = electron.get_number_converged()
print "number of states", last_e

plot(PE, title="Potential energy")
for i in range(0,3):
    plot_wave(electron, i, Q, title="electron "+str(i))

plot_wave(electron, last_e-1, Q,title="electron "+ str(last_e-1))
interactive()


-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.