Quantum Variational Solver

This small script allows you to calculate the approximate ground state energy of a Hamiltonian given a wavefunctional form inputted as a symbolic expression from sympy. The energy minimization method is not the best, and especially questionable with more than one parameter, but it works as a toy model to demonstrate how variational methods are used in QM. There are several common forms defined in the Forms.py file, which look something like:

# Spherical Exponential (radial but normalized over theta and phi)
def SphericalExponential():
    # (1)
    r, b, theta, phi = sp.symbols(r'r b \theta \phi', real=True, positive=True)

    # (2)
    var = [r, theta, phi]
    bounds = [[0, sp.oo], [0, sp.pi], [0, 2*sp.pi]]
    volumeElement = r**2 * sp.sin(theta)
    params = b

    # (3)
    form = sp.exp(-b * r)

    # (4)
    return form, var, bounds, volumeElement, params

1. First, we declare some some symbols for both variables and parameters. It is important to be as specific with these as possible, which is why I have specified that they are both real and positive. Without these constraints, you may run into some issues parameterizing the hamiltonian/ntegrating.

2. Next up, we make lists of the variables that will be integrated over, the integral bounds for normalizing/evaluating expectation values, the volume element for integration, and the parameters that can be varied.

3. This defines the actual form the wavefunction, and can include any of the parameters or variables from (1)

4. We return all of the information so that it can be used to parameterize the hamiltonian and evaluate expectation values.

Note that you do not need to explicitly normalize a form, as it can be done with

form, var, bounds, volumeElement, params = SphericalExponential() # Or any other form

normForm = Variation.normalize(form, var, bounds, volumeElement)

After this, we parameterize the hamiltonian in terms of any of the variational parameters in params:

energyFunc = Variation.parameterizeHamiltonian(Hamiltonian, normForm, var, bounds, volumeElement)

Some example Hamiltonians are defined in the Jupyter Notebook GenericSolver.ipynb, but they should look something like:

def BohrHamiltonian(psi, var):
    r, theta, phi = var
    laplacian = 1/r**2 * sp.diff(r**2 * sp.diff(psi, r), r)
    kineticTerm = -hbar**2 / (2*me) * laplacian
    potentialTerm = -e**2 / (4 * e0 * sp.pi * r) * psi
    return (kineticTerm + potentialTerm)/e

We unpack the variables contained in var (which should be a tuple for forms with more than one variable), and then do whatever is defined by the hamiltonian. Since operators may be derivatives, we have to define the hamiltonian as acting on a wavefunction psi, though it may not always be changed (eg. the potential term above just multiples by it since 1/r doesn't act on the wavefcuntion).

Finally, we pass the parameterized energy function as well as what we are optimizing with respect to, and hopefully we'll find a minimum in the curve.

optimalParam, minEnergy = Variation.minimizeHamiltonian(energyFunc, params, plot=True)

The first return is the (possibly a list) of optimal parameter values, and the second is the energy function evaluated at those points. This function mostly uses sympy's solve function to find the zero of the derivative of the function, but this only works sometimes. I have implemented a basic multidimensional optimizer which iterates over each parameter and makes several passes, but it doesn't really work all that well.

Results

If you'd like to compare the symbolic solver to an analytical derivation of the ground state energy of the hydrogen atom (with relativistic corrections) using a radial exponential, see GroundStateHydrogen.pdf (derivation) and GroundStateHydrogen.ipynb (plotting). The two methods line up very nicely, and give results that are very close to the true ground state energy -- within the error generated by having only so many sig figs for physical constants. This is especially a problem since some of the constants, like the Planck constant, are raised to the 4th power.

Calculated 1/b: 5.29161334475323e-11 m
Analytic Bohr radius: 5.29177210903e-11 m
Variational relativistic correction: -0.0001811521353554502 eV
Analytic relativistic correction: -0.0009056519849917295 eV

References

Griffiths, D. J., & Schroeter, D. F. (2018). Introduction to quantum mechanics (Third edition). Cambridge University Press.

Meurer A, Smith CP, Paprocki M, Čertík O, Kirpichev SB, Rocklin M, Kumar A, Ivanov S, Moore JK, Singh S, Rathnayake T, Vig S, Granger BE, Muller RP, Bonazzi F, Gupta H, Vats S, Johansson F, Pedregosa F, Curry MJ, Terrel AR, Roučka Š, Saboo A, Fernando I, Kulal S, Cimrman R, Scopatz A. (2017) SymPy: symbolic computing in Python. PeerJ Computer Science 3:e103 https://doi.org/10.7717/peerj-cs.103