Combining two potentials in python

Question

So I need to define the potential in python. I have infinite square well potential, and it has, in the middle the potential barrier.

For the barrier I have this:

    V_x = zeros([Npts],float)
for i in range(0,Npts,1):
    if x[i] > 0 and x[i]<width:
        V_x[i]=v0

Npts is the the length of the x, and x is defined as array of increasing values from xmin to xmax (x=arange(Xmin,Xmax+0.001,dx)).

How to include the infinite potential into this?

It has some length (let's say from -100 to 100) and it has to behave like impenetrable wall (the function at the -100 and 100 has to be zero).

I could combine the barrier with this potential separately (e.g. V_b for the barrier, and V_p for the infinite potential and then the V_x could be V_x=V_b+V_p).

Do you think that could work? I have the problem defining that infinite potential...

Answer 1

I would suggest just using a number for your infinite potential that is a few orders (3-4) of magnitude larger than your finite potential. It will not make much of a difference in the energies of your states to make it any larger. In fact, you'll quickly be limited more by the precision of the floats used when solving the Hamiltonian. You really only need to set the ends to be 'infinite' so something like the following should work: V_x[0] = V_x[-1] = abs(v0)*10**4. Of course, you can set it to be larger if you want, but I don't think that it will make much of a difference. If you aren't setting the ends to be infinite, then why are you trying to solve for the value of your wave function in large portions of the region in which you know that the wave function will be zero?

Also, be careful that your infinite potential isn't of compararable magnitude to h-bar^2/2/m/delta_x^2 in the units that you are using. It's all making the 'infinite' potential be a few orders of magnitude larger than both the finite potential and the kinetic energy terms in your Hamiltonian.

Answer 2

Numpy has a Inf float variable, so depending on how you are using the potential, you can define it explicitly:

import numpy as np
Npts = 100
v0 = 10.0
V_x = v0*np.ones((Npts,),dtype=float)
x = np.linspace(-20,20,Npts)
width = 15.0
# now replace those values for the potential in the impenetrable
# region with np.Inf 
ii = ((x>width) | (x<0)).nonzero()
V_x[ii] = np.Inf

This gives you a system with x between -20 and 20, where V_x = v0 for 0 < x < 15 and V_x = Inf for x > 15 or x < 0

Note: As a general suggestion, there are also some more efficient ways to set up numpy arrays rather than using for loops, as is shown in the above code.