Calculating eigenstates and -energies of operators

This web page can be downloaded as notebook: eigenstates.ipynb (Jupyter) or eigenstates.md (Markdown)

There are a number of use cases where we want to calculate the eigenstates of a Hamiltonian, for example, if we start a calculation from an ensemble of low-temperature states. Here, we show some basic and advanced use cases around energy eigenstates.

Note that there exist other approaches for calculating eigenstates of the Hamiltonian. In particular, for few low-energy eigenstates, you can try imaginary-time propagation, as detailed in Ground-state relaxation and finite temperatures.

Basic eigenstate calculation

As example system, we choose a modified double-well potential (Razavy potential). Diagonalizing the Hamiltonian is done with wp.diagonalize(), which yields the eigenstates and -energies sorted by the latter.

import numpy as np
import wavepacket as wp

def razavy_potential(dvr_points):
    val = np.cosh(dvr_points)
    return -0.7 * val + 0.01 * val**2
    
grid = wp.grid.Grid(wp.grid.PlaneWaveDof(-7, 7, 128))
hamiltonian = wp.operator.CartesianKineticEnergy(grid, 0, 0.5) + wp.operator.Potential1D(grid, 0, razavy_potential)

for energy, state in wp.diagonalize(hamiltonian):
    # placeholder for doing something with the eigenstate
    print(f"E = {energy:.4} a.u., |psi|^2 = {wp.trace(state):.4}")

E = -9.002 a.u., |psi|^2 = 1.0
E = -9.002 a.u., |psi|^2 = 1.0
E = -4.012 a.u., |psi|^2 = 1.0
E = -4.011 a.u., |psi|^2 = 1.0
E = -1.16 a.u., |psi|^2 = 1.0
E = -0.9569 a.u., |psi|^2 = 1.0
E = 0.2135 a.u., |psi|^2 = 1.0
E = 1.374 a.u., |psi|^2 = 1.0
E = 2.813 a.u., |psi|^2 = 1.0
E = 4.463 a.u., |psi|^2 = 1.0
E = 6.305 a.u., |psi|^2 = 1.0
E = 8.324 a.u., |psi|^2 = 1.0
E = 10.51 a.u., |psi|^2 = 1.0
E = 12.86 a.u., |psi|^2 = 1.0
E = 15.37 a.u., |psi|^2 = 1.0
E = 18.03 a.u., |psi|^2 = 1.0
E = 20.84 a.u., |psi|^2 = 1.0
E = 23.79 a.u., |psi|^2 = 1.0
E = 26.88 a.u., |psi|^2 = 1.0
E = 30.11 a.u., |psi|^2 = 1.0
E = 33.48 a.u., |psi|^2 = 1.0
E = 36.99 a.u., |psi|^2 = 1.0
E = 40.62 a.u., |psi|^2 = 1.0
E = 44.39 a.u., |psi|^2 = 1.0
E = 48.28 a.u., |psi|^2 = 1.0
E = 52.3 a.u., |psi|^2 = 1.0
E = 56.45 a.u., |psi|^2 = 1.0
E = 60.72 a.u., |psi|^2 = 1.0
E = 65.11 a.u., |psi|^2 = 1.0
E = 69.62 a.u., |psi|^2 = 1.0
E = 74.26 a.u., |psi|^2 = 1.0
E = 79.01 a.u., |psi|^2 = 1.0
E = 83.88 a.u., |psi|^2 = 1.0
E = 88.87 a.u., |psi|^2 = 1.0
E = 93.97 a.u., |psi|^2 = 1.0
E = 99.19 a.u., |psi|^2 = 1.0
E = 104.5 a.u., |psi|^2 = 1.0
E = 110.0 a.u., |psi|^2 = 1.0
E = 115.5 a.u., |psi|^2 = 1.0
E = 121.2 a.u., |psi|^2 = 1.0
E = 127.0 a.u., |psi|^2 = 1.0
E = 132.9 a.u., |psi|^2 = 1.0
E = 138.9 a.u., |psi|^2 = 1.0
E = 145.0 a.u., |psi|^2 = 1.0
E = 151.2 a.u., |psi|^2 = 1.0
E = 157.5 a.u., |psi|^2 = 1.0
E = 164.0 a.u., |psi|^2 = 1.0
E = 170.5 a.u., |psi|^2 = 1.0
E = 177.1 a.u., |psi|^2 = 1.0
E = 183.9 a.u., |psi|^2 = 1.0
E = 190.7 a.u., |psi|^2 = 1.0
E = 197.7 a.u., |psi|^2 = 1.0
E = 204.8 a.u., |psi|^2 = 1.0
E = 211.9 a.u., |psi|^2 = 1.0
E = 219.2 a.u., |psi|^2 = 1.0
E = 226.6 a.u., |psi|^2 = 1.0
E = 234.0 a.u., |psi|^2 = 1.0
E = 241.6 a.u., |psi|^2 = 1.0
E = 249.3 a.u., |psi|^2 = 1.0
E = 257.0 a.u., |psi|^2 = 1.0
E = 264.9 a.u., |psi|^2 = 1.0
E = 272.9 a.u., |psi|^2 = 1.0
E = 280.9 a.u., |psi|^2 = 1.0
E = 289.1 a.u., |psi|^2 = 1.0
E = 297.4 a.u., |psi|^2 = 1.0
E = 305.7 a.u., |psi|^2 = 1.0
E = 314.2 a.u., |psi|^2 = 1.0
E = 322.7 a.u., |psi|^2 = 1.0
E = 331.4 a.u., |psi|^2 = 1.0
E = 340.2 a.u., |psi|^2 = 1.0
E = 349.0 a.u., |psi|^2 = 1.0
E = 357.9 a.u., |psi|^2 = 1.0
E = 367.0 a.u., |psi|^2 = 1.0
E = 376.1 a.u., |psi|^2 = 1.0
E = 385.3 a.u., |psi|^2 = 1.0
E = 394.6 a.u., |psi|^2 = 1.0
E = 404.1 a.u., |psi|^2 = 1.0
E = 413.6 a.u., |psi|^2 = 1.0
E = 423.2 a.u., |psi|^2 = 1.0
E = 432.9 a.u., |psi|^2 = 1.0
E = 442.7 a.u., |psi|^2 = 1.0
E = 452.5 a.u., |psi|^2 = 1.0
E = 462.5 a.u., |psi|^2 = 1.0
E = 472.6 a.u., |psi|^2 = 1.0
E = 482.7 a.u., |psi|^2 = 1.0
E = 493.0 a.u., |psi|^2 = 1.0
E = 503.3 a.u., |psi|^2 = 1.0
E = 513.7 a.u., |psi|^2 = 1.0
E = 524.3 a.u., |psi|^2 = 1.0
E = 534.9 a.u., |psi|^2 = 1.0
E = 545.6 a.u., |psi|^2 = 1.0
E = 556.4 a.u., |psi|^2 = 1.0
E = 567.2 a.u., |psi|^2 = 1.0
E = 578.2 a.u., |psi|^2 = 1.0
E = 589.3 a.u., |psi|^2 = 1.0
E = 600.4 a.u., |psi|^2 = 1.0
E = 611.7 a.u., |psi|^2 = 1.0
E = 623.0 a.u., |psi|^2 = 1.0
E = 634.4 a.u., |psi|^2 = 1.0
E = 645.9 a.u., |psi|^2 = 1.0
E = 657.5 a.u., |psi|^2 = 1.0
E = 669.1 a.u., |psi|^2 = 1.0
E = 680.9 a.u., |psi|^2 = 1.0
E = 692.7 a.u., |psi|^2 = 1.0
E = 704.8 a.u., |psi|^2 = 1.0
E = 716.9 a.u., |psi|^2 = 1.0
E = 729.2 a.u., |psi|^2 = 1.0
E = 741.2 a.u., |psi|^2 = 1.0
E = 752.8 a.u., |psi|^2 = 1.0
E = 763.4 a.u., |psi|^2 = 1.0
E = 776.1 a.u., |psi|^2 = 1.0
E = 789.4 a.u., |psi|^2 = 1.0
E = 810.4 a.u., |psi|^2 = 1.0
E = 827.1 a.u., |psi|^2 = 1.0
E = 831.6 a.u., |psi|^2 = 1.0
E = 940.8 a.u., |psi|^2 = 1.0
E = 943.3 a.u., |psi|^2 = 1.0
E = 1.101e+03 a.u., |psi|^2 = 1.0
E = 1.103e+03 a.u., |psi|^2 = 1.0
E = 1.311e+03 a.u., |psi|^2 = 1.0
E = 1.314e+03 a.u., |psi|^2 = 1.0
E = 1.58e+03 a.u., |psi|^2 = 1.0
E = 1.585e+03 a.u., |psi|^2 = 1.0
E = 1.922e+03 a.u., |psi|^2 = 1.0
E = 1.931e+03 a.u., |psi|^2 = 1.0
E = 2.361e+03 a.u., |psi|^2 = 1.0
E = 2.364e+03 a.u., |psi|^2 = 1.0
E = 3.009e+03 a.u., |psi|^2 = 1.0

Note that the direct diagonalization does not do fancy tricks. It only extracts the (dense) matrix representation of the operator, diagonalizes it using numpy.linalg.eigh, and transforms the output into Wavepacket data structures. It is therefore rather robust, but not terribly efficient.

Computation times of such standard solvers scale with the third power of the number of grid points, memory with the square of the number of grid points. At a size of 16 bytes per complex value, you need maybe 1 GB of memory for 10,000 grid points (100x100 grid), also for work storage, and several minutes of computation time. This is about the maximum size that you can reasonably use this procedure for.

Advanced usage: Largest and smallest eigenstates

In some applications, you are not interested in all eigenstates, but only some of them. Usually, you start by wrapping the calculation results in a list. This is slightly inefficient if you do not need all states, but the wasted effort is much smaller than the effort for the eigenvalue calculation in the first place.

results = list(wp.diagonalize(hamiltonian))
len(results)

We now have a list of (energy, state) pairs, and can manipulate them with usual Python list comprehensions. For example, to get all energies up to a cut-off,

results_with_negative_energy = [(energy, state) for energy, state in results if energy < 0]
len(results_with_negative_energy)

The same approach can be used with slicing, sorting etc. The resulting list can always be iterated over like the initial generator

sliced_results = results[0:10]

inverted_order = sliced_results
inverted_order.sort(reverse=True)

for energy, state in inverted_order:
    print(f"E = {energy:.4} a.u., |psi|^2 = {wp.trace(state):.4}")

E = 4.463 a.u., |psi|^2 = 1.0
E = 2.813 a.u., |psi|^2 = 1.0
E = 1.374 a.u., |psi|^2 = 1.0
E = 0.2135 a.u., |psi|^2 = 1.0
E = -0.9569 a.u., |psi|^2 = 1.0
E = -1.16 a.u., |psi|^2 = 1.0
E = -4.011 a.u., |psi|^2 = 1.0
E = -4.012 a.u., |psi|^2 = 1.0
E = -9.002 a.u., |psi|^2 = 1.0
E = -9.002 a.u., |psi|^2 = 1.0