quantum-visualizations/qmsolve

⚛️ A module for solving and visualizing the Schrödinger equation.

qmsolve's Issues

Trying To Specify a 3D Potential, Wierd Gaps

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 :)

Simulate make Ising Hamiltonin

Helo,
I try make Ising Hamiltonian simulation, how represent H for simulation with two qubits?

Thanks you

Request: Add Eigsh-Cupy Mode

CuPy also has eigsh and it makes computation much faster, the code should be almost identical to what is alraedy there

        if method == 'eigsh':
            from scipy.sparse.linalg import eigsh

            # Note: uses shift-invert trick for stability finding low-lying states
            # Ref: https://docs.scipy.org/doc/scipy/reference/tutorial/arpack.html#shift-invert-mode

            eigenvalues, eigenvectors = eigsh(H, k=max_states, which='LM', sigma=min(0, self.E_min))

Simply change the import (and condition) and it should work.

Solving the Schrodinger Equation in SI untis.

It is not clear from the examples how to solve the Schrodinger equation in SI units from the examples. Could you please clarify on this?

Bohmian trajectories for time dependent wave functions

I've been studying the 2D time-dependent two-slit diffraction example. I'm trying to understand the units of the time degree of freedom. If I modify the potential to be all zero, then I should just get the free particle Gaussian wave packet in 2D which is solvable analytically. When I do this, it looks like the phase velocity at the center of the Gaussian wave packet is quite a bit slower than the motion of the wave packet's density center, and I think they should be about the same. I'm trying to use qmsolve to simulate Bohmian trajectories and eventually stochastic mechanics trajectories, and so I want to make sure that the Gaussian case looks right as a sanity check. Peter Holland's book "The Quantum Theory of Motion", chapter 4, covers this in some detail. Has anybody thought to try and use qmsolve for calculating Bohmian trajectories yet? It seems like an obvious thing to do.

Non-printable characters

Some code contains non-printable characters that causes problems when trying to run simulations.
Manually solving this manually is tedious. Any chance of using only printable characters?

Animation error

The simulation seems to run ok. However a plot at t=0 is plotted but no animation.

The following warning is displayed:

UserWarning: Animation was deleted without rendering anything. This is most likely not intended. To prevent deletion, assign the Animation to a variable, e.g. anim, that exists until you have outputted the Animation using plt.show() or anim.save().
warnings.warn(

I have tried using plt.show() or anim.save(). but no luck.

The solver doesn't work correctly when there are negative eigenvalues

Code to reproduce:

import numpy as np
from qmsolve import Hamiltonian, SingleParticle, init_visualization


#interaction potential
def single_gaussian_well(particle):
	𝜇 = 0.0
	σ = 0.5
	V = -600* np.exp((-(particle.x)**2 -(particle.y-𝜇)**2 ) / (2*σ**2))
	return V



H = Hamiltonian(particles = SingleParticle(), 
				potential = single_gaussian_well, 
				spatial_ndim = 2, N = 100, extent = 8)


eigenstates = H.solve(max_states = 60)

print(eigenstates.energies)
visualization = init_visualization(eigenstates)
#visualization.plot_eigenstate(6)
visualization.slider_plot()

The problem is that :
eigsh(H, k = max_states, which='LM', sigma=0)
doesn't work with negative eigenvalues.

Because
eigsh(H, k = max_states, which='SA')
works with positive eigenvalues, I suggest adding an option to choose the correct solver, or just scaling the potential to be always positive.

Energy Visualization problem in 4 Gaussian Wells Example

I am having some issues learning how to work with this package. Some of the examples don't seem to work the same way on my computer.
Most work, but there are two that don't.

Specifically, 3D_two_gaussian_wells.py, line 21, references a variable "N0" which doesn't exist, leading to an error.

Also, I am having problems with the animation in 4 Gaussian Wells. I see the wavefunction visualization, but I don't see the energy levels? it's odd, because it's otherwise exactly the same as the example.

Request for Final year project

How are you? I have seen your videos on YouTube. Physics simulation videos are very beautiful and I am a final semester student of bachelor degree in physics and my final year project is also on physics simulation. It didn't work. Can you help me? Make me project at reasonable price.I share more details of the project with you if you want to make it.

Imaginary part of the wave function

In the time independent harmonic oscillator from example, the eigen states doesn't have any imaginary part contrary to what could be expected by the analytical resolution of the harmonic oscillator. (exp(imphi))

Do you have any idea why, and how I could fix it ?

Visualization of potential isn't correct for negative values of V0 in 1D_potential_barrier.py

I am running the example 1D_potential_barrier.py. In the function potential barrier, if set V0 to a negative value then the visualization of the potential isn't correct. For negative values of V0 the potential is -1 when it should be 0 and -2 when it should be the V0.

I think the issue lies in line 432 of the animate function.

In 1D is time double faster than p_x0 in initial gaussians?

Hello, when running the sample 1d potential barrier with 0 pot value, if I give 50 AU v0 speed(same that p_x0 because mass is 1), it takes 0.55 AU to walk 50 AU, when it have to take 1 AU.
Can this be fixed replacing?:
np.exp(p_x0*particle.x*1j)
with
np.exp(p_x0*particle.x*0.5j)
Take note that I've modified your code for plotting AU instead of femtoseconds:
time_ax.set_text(u"t = {} au".format("%.3f" % ((animation_data['t'] * TAUFMS))))
and the same for space Amstrong plot