This branch implements an alternative eigensolver routine that has much quicker convergence in some circumstances. The new solver is based on scipy.sparse.linalg.lobpcg. By default the code behaves as before and uses the eigsh solver. To use the new method pass method="lobpcg" to Hamiltonian.solve.

Unfortunately, lobpcg is not nearly as reliable as eigsh. In some cases, lobpcg is very slow and produces very wrong results. I have seen better convergence in situations that do not require large grid size. I have only tested on harmonic oscillator. I see no benefit for 1 or 2 dimensions, and a large benefit in 3 dimensions. I have written a script to demonstrate this as discussed below. Further tests are needed to determine the stability in a wider range situations.

Though not strictly required, the implementation benefits from providing a preconditioning operator that approximates the inverse of the matrix to be diagonalized (the hamiltonian in this case). The actual inverse cannot be efficiently computed. The current code uses the "Jacobi" or "diagonal" conditioning matrix which is a matrix of zeros except along the diagonal. On the diagonal, the values are simply 1/diagonal of the original hamiltonian values.

To demonstrate the usage, I have created a new example script in a new examples/ folder (btw, perhaps we should put all examples here; top level is looking very crowded :). The example script solves for the eigenstates of a single particle in a 3d isotropic HO. The lobpcg solver is ~40 times faster in my test with relative errors <2% for the first 5 eigenvalues. The 10th eigenvalue has 25%. Here is an example output from running the script:

------------------------------
Solving with scipy.sparse.linalg.eigsh
Computing...
Took 36.023282051086426
Eigenvalues (eigsh) : [ 6.03533971 10.04416034 10.04416034 10.04416034 14.0257434  14.0257434
 14.0257434  14.05298098 14.05298098 14.05298098]
------------------------------
Solving with scipy.sparse.linalg.lobpcg
Computing...
Took 0.9391400814056396
Eigenvalues (lobpcg): [ 6.03600915 10.05316063 10.0997381  10.24087437 14.05243788 14.2446791
 14.31203711 14.66924361 15.23159918 18.16693888]
------------------------------
Solving with foo
Computing...
foo solver has not been implemented. Use one of ('eigsh', 'lobpcg')
------------------------------
Diff. over mean: [1.10914213e-04 8.95670917e-04 5.51807386e-03 1.93949903e-02
 1.90143947e-03 1.54886751e-02 2.02057964e-02 4.29119016e-02
 8.04941167e-02 2.55367358e-01]

Process finished with exit code 0

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.

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

Include the superpositions functions from the main branch in the prototype as well. This requires adding the new class ComplexSliderWidget. I've modified 1D_harmonic_oscillator.py and 2D_harmonic_oscillator.py to present examples of the new method.

Following guidance here: https://docs.scipy.org/doc/scipy/reference/tutorial/arpack.html#shift-invert-mode. Simple tests show dramatic speedup in eigensolver.

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.

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

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?

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.