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

[rafael-fuente]

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.

[dhudsmith]

One easy fix would be to always set sigma to the minimum value of the potential evaluated on the grid. Using the shift-invert trick with which="LM", the solver will find eigenvalues close to sigma.

[rafael-fuente]

Fixed