Trying To Specify a 3D Potential, Wierd Gaps

[cgbsu]

Hello,

Your library is amazing and so are your examples. I was working on something similar for my research then stumbled upon your project.

I am trying to specify a simple barrier for Quantum Tunneling. I copy/pasted the method for rendering eigen modes and used it to render the potential

I specify the potential with this code

def tunnelingBarrier(
            particle, potential = 1, 
            offsetX = .2, offsetY = .5, offsetZ = .5, 
            width = .2, height = 1, depth = 1
        ):
    extent = np.abs(particle.x.max()) + np.abs(particle.x.min())
    offsetX *= particle.x.max()
    width *= particle.x.max()
    return np.where(
            (particle.x < (offsetX + width)) & (particle.x > offsetX), 
            np.ones(particle.x.shape) * potential, 
            np.zeros(particle.x.shape), 
        )

Here is what I get:


It took quite a lot of struggle getting too this point, and thankfully it is sort of working
I can just see a an evanescent probability density on the other side of where the barrier might be.

However its quite wonky, first there is the matter of the gap, the value should be constant wherever the conditions are satisfied.

If I double the offset (without changing particle.x.max()) the width of the barrier increases as well, while it should not be effected. Its like there are 2 of them.

I did modify the existing way of rendering things, so maybe I did something wrong there and its just a graphical thing?

from mayavi import mlab
import numpy as np
from qmsolve.util.constants import *

def plot_potential(hamiltonian, contrast_vals= [0.1, 0.25], plot_type = "volume"):
    potential = hamiltonian.Vgrid
    mlab.figure(1, bgcolor=(0, 0, 0), size=(700, 700))

    abs_max= np.amax(np.abs(potential))
    potential = (potential)/(abs_max)

    L = hamiltonian.extent / 2 / Å
    N = hamiltonian.N

    vol = mlab.pipeline.volume(mlab.pipeline.scalar_field(potential))

    # Change the color transfer function
    from tvtk.util import ctf
    c = ctf.save_ctfs(vol._volume_property)
    c['rgb'] = [[-0.45, 0.3, 0.3, 1.0],
                [-0.4, 0.1, 0.1, 1.0],
                [-0.3, 0.0, 0.0, 1.0],
                [-0.2, 0.0, 0.0, 1.0],
                [-0.001, 0.0, 0.0, 1.0],
                [0.0, 0.0, 0.0, 0.0],
                [0.001, 1.0, 0.0, 0.],
                [0.2, 1.0, 0.0, 0.0],
                [0.3, 1.0, 0.0, 0.0],
                [0.4, 1.0, 0.1, 0.1],
                [0.45, 1.0, 0.3, 0.3]]

    c['alpha'] = [[-0.5, 1.0],
                  [-contrast_vals[1], 1.0],
                  [-contrast_vals[0], 0.0],
                  [0, 0.0],
                  [contrast_vals[0], 0.0],
                  [contrast_vals[1], 1.0],
                 [0.5, 1.0]]
    if plot_type == 'volume':
        ctf.load_ctfs(c, vol._volume_property)
        # Update the shadow LUT of the volume module.
        vol.update_ctf = True

        mlab.outline()
        mlab.axes(xlabel='x [Å]', ylabel='y [Å]', zlabel='z [Å]',nb_labels=6 , ranges = (-L,L,-L,L,-L,L) )
        #azimuth angle
        φ = 30
        mlab.view(azimuth= φ,  distance=N*3.5)
        mlab.show()

    if plot_type == 'abs-volume':
    
        abs_max= np.amax(np.abs(potential))
        psi = (potential)/(abs_max)

        L = hamiltonian.extent/2/Å
        N = hamiltonian.N

        vol = mlab.pipeline.volume(mlab.pipeline.scalar_field(np.abs(psi)), vmin= contrast_vals[0], vmax= contrast_vals[1])
        # Change the color transfer function

        mlab.outline()
        mlab.axes(xlabel='x [Å]', ylabel='y [Å]', zlabel='z [Å]',nb_labels=6 , ranges = (-L,L,-L,L,-L,L) )
        #azimuth angle
        φ = 30
        mlab.view(azimuth= φ,  distance=N*3.5)
        mlab.show()

    elif plot_type == 'contour':
        psi = potential
        L = hamiltonian.extent/2/Å
        N = hamiltonian.N
        isovalue = np.mean(contrast_vals)
        abs_max= np.amax(np.abs(potential))
        psi = (psi)/(abs_max)

        field = mlab.pipeline.scalar_field(np.abs(psi))

        arr = mlab.screenshot(antialiased = False)

        mlab.outline()
        mlab.axes(xlabel='x [Å]', ylabel='y [Å]', zlabel='z [Å]',nb_labels=6 , ranges = (-L,L,-L,L,-L,L) )
        colour_data = np.angle(psi.T.ravel())%(2*np.pi)
        field.image_data.point_data.add_array(colour_data)
        field.image_data.point_data.get_array(1).name = 'phase'
        field.update()
        field2 = mlab.pipeline.set_active_attribute(field, 
                                                    point_scalars='scalar')
        contour = mlab.pipeline.contour(field2)
        contour.filter.contours= [isovalue,]
        contour2 = mlab.pipeline.set_active_attribute(contour, 
                                                    point_scalars='phase')
        s = mlab.pipeline.surface(contour, colormap='hsv', vmin= 0.0 ,vmax= 2.*np.pi)

        s.scene.light_manager.light_mode = 'vtk'
        s.actor.property.interpolation = 'phong'


        #azimuth angle
        φ = 30
        mlab.view(azimuth= φ,  distance=N*3.5)

        mlab.show()

The field is normalized but if its constant over a region then it should have the same color/value throughout the region. I looked at the code for Hamiltonian, the generation of the potential looks pretty straightforward, I'm not sure what is causing it.

So I would like to ask:
A: How do I resolve this, and is it me or qmsolve?
B: Could we have some documentation on how to specify potentials?

Thank you :)

[cgbsu]

P.s I can confirm they are changing size + position. Maybe there is some black parts of the thin barrier on the outsides too?