Solving the Helium Eigenvalue Problem—where to start?

Specific Domains
, , ,
1

I have inherited Fortran code (the same one in this paper) that solves the time-independent Schrodinger equation for the ground state of the helium atom. This is generally a very difficult problem because the Schrodinger equation for helium is a non-separable 6-D PDE:

\hat{H}\psi(\vec{r_1}, \vec{r_2}) = E \psi(\vec{r_1}, \vec{r_2}),

where E is a constant corresponding to the total energy of the system, and \hat{H} is an operator defined as follows:

\nabla^2_1 + \nabla^2_2 - \frac{2}{r_1} - \frac{2}{r_2} + \frac{1}{|\vec{r_1} - \vec{r_2}|}

The Fortran code that I have reduces this to a 2-D problem by expanding \psi in terms of coupled spherical harmonics, then summing over the different parameters of those harmonics—this is referred to as the method of partial waves. The two radial coordinates r_1 and r_2 are then mapped to a variable finite-element grid. Finally, the program solves the problem by solving the time-dependent Schrodinger equation (TDSE) instead (\hat{H}\psi(\vec{r_1}, \vec{r_2}) = i \frac{\partial}{\partial t} \psi(\vec{r_1}, \vec{r_2})), but with an imaginary time parameter. This is because if you naively solve the TDSE by separating the derivative, you obtain the following:

\psi(\vec{r_1}, \vec{r_2}, t + dt) = e^{i \hat{H}dt} \psi(\vec{r_1}, \vec{r_2}, t)

If you change t to \tau = it, then the matrix exponential above becomes real and decaying. If \psi is a convenient test function (say, a Gaussian) that contains a mix of all possible eigenfunctions, then since e^{i \hat{H}d\tau} returns e^{i E_n d\tau} when applied to any given eigenfunction \psi_n, those states with higher E_n decay to zero more quickly as the wavefunction is propagated, leaving only the lowest-E_n state, which is the ground state. In-code, this propagation is done via the split-operator method, Fourier-transforming the system such that each operator matrix is diagonal when applied, except for the coupling term \frac{1}{|\vec{r_1} - \vec{r_2}|} . The code I have currently takes a couple minutes to run on 100 cores on an HPC cluster and returns both the right ground state energy (-2.903 au) and the ground-state wavefunction.

So, if I have code that works, runs quickly, and I know how to use it, why am I trying to reinvent the wheel? I’m glad you asked, imaginary curious interlocutor. There’s several reasons:

  1. The code is nigh-unmaintainable. There’s minimal comments, obscure-at-best variable names, and this is on its 4th generation of being handed down, so I can’t easily ask questions about how it works or what effects making changes would have.
  2. The code is 20 years old. This is no bad thing in itself, but in the intervening decades, as computers have gotten better and people have put more and more time into efficiently solving difficult PDEs, I’d like to explore what new methods are available to me.
  3. I want to understand how my own code works. This is more about my own learning than practicality, but I’d like to move beyond plugging-and-chugging someone else’s code and towards understanding how to actually solve these kinds of problems, so that I can better extend it and explain it.

I know that Julia has many options for solving PDEs, but it’s hard for me to figure out which methods are right for solving this kind of problem, especially since the choices made for the code I have were determined by constraints 20 years ago. Finite differences were too inaccurate and spectral methods took too long, but I don’t know if that’s changed. For example, MethodOfLines.jl says it can handle spherical Laplacians, but I can’t find further documentation on that. None of the packages I’ve looked at (MethodOfLines, Ferrite, Gridap) talk about coupled terms like the one I’m dealing with. I would greatly appreciate some direction from one of the fine linear algebra experts on this forum. Thank you!

2

In numerics this would typically be called a “spectral” method. (It’s not really “2d” in the sense that the number of spherical harmonics that you need is still determined by the 3d resolution, but you can get exponentially fast convergence with the number of harmonics.)

My understanding is that “imaginary time evolution” like this has mostly been superseded by implicitly restarted Lanczos, LOBPCG, and similar Krylov-type methods.

You can still use this in a Krylov method in order to do your matrix–vector multiplications. There are packages like FastSphericalHarmonics.jl which can help.

The method of “partial waves” is a spectral method.

3

Thank you for the suggestions! With the current code, 6 partial waves is enough for convergence to the ground state, so that tracks with what you’re saying.

Also, my fault for not proofreading/checking my sources thoroughly enough—I meant to say global spectral methods were considered too computationally demanding relative to the level of resolution needed, hence the use of finite elements (in this case, spanned by Legendre polynomials, so the grid points correspond to the roots of the polynomials for easy quadrature). I do think part of the difficulty is that I never had a proper class on linear algebra or PDE solving, so I’m very much picking things up as I go.

4

That sounds like a pseudo-spectral collocation method? i.e. you use a spectral basis for the unknowns, and you enforce the equations at a grid of sample points, carefully chosen (usually related roots of the spectral basis or its derivatives).

5

That sounds right, though I haven’t encountered that exact phrasing before.

6

Some years ago I did something here:

which is the Julia version of the fortran code used in

maybe there is something useful there, if not, apologies for the noise.

Edit: There is not a a lot of documentation or even a proper README, but it works, maybe I should bring this back to life :D.

7

So, thinking about the problem a bit more and doing some more research, I’ve determined the following:

  • Yes, the method used in the paper is a pseudospectral collocation method (per Wikipedia, apparently pseudospectral is a synonym for discrete variable representation, which is indeed what the Fortran code uses. Seems I’m one of today’s lucky ten thousand.)
  • I inadvertently left a detail out of my description; after expanding \psi and 1/|\vec{r_1} - \vec{r_2}| in terms of coupled spherical harmonics, the entire equation is integrated over both pairs of angular coordinates, which (because of the orthogonality properties of the coupled harmonics) yields the following, where V_{F0}(i,j) is the expanded and integrated form of 1/|\vec{r_1} - \vec{r_2}|, and i and j index bundles of the integer partial wave parameters L, l_1, and l_2:
-\frac{1}{2}\frac{\partial^2}{\partial r_1^2} -\frac{1}{2}\frac{\partial^2}{\partial r_2^2} + \frac{l_1(l_1 + 1)^2}{2r_1^2} + \frac{l_2(l_2 + 1)^2}{2r_2^2} - \frac{2}{r_1} - \frac{2}{r_2} \\+ \sum_j V_{F0}(i,j) = i \frac{\partial}{\partial t} \psi_i(r_1, r_2)
  • The above coupled PDEs are what are properly discretized on a 2D grid of Legendre polynomial roots and solved via reverse imaginary time propagation.
  • If I’m purely interested in the ground-state eigenvalues and eigenvectors of this problem, then if I can construct the operators as a matrix, I can use KrylovKit.jl or ArnoldiMethod.jl to solve the eigenvalue problem for a given set of partial waves, making sure that I use enough partial waves to converge on the correct ground-state energy.

So, with this in mind, I think I can ask some more specific questions, namely:

  1. Is pseudospectral collocation (what the paper calls finite element-discrete variable representation) the only viable approach, or is it worth trying a different method?
  2. Based on the answer to the above, what Julia packages will help me specify the grid and construct the operators for this problem? I’m not really sure what the operators will look like if I’m not doing the split-operator + reverse imaginary time method.
8

Took me a minute to think through but I think the split-operator method just drops out fully.

split-operator method

Split-operator method approximates the time-evolution operator if the Hamiltonian can be separated into 2 non-commuting parts H=H_A + H_B that are simple on their own:

U(t) = exp(-i H t) = exp(-iH \Delta t)^N\approx [exp(-iH_A \Delta t) exp(-iH_B \Delta t)]^N

In code this looks something like:

function U(t, psi, N)
    diagonal_H_A = ...
    diagonal_H_B = ...

    diagonal_U_A = @. exp(-im*diagonal_H_A*t/N)
    diagonal_U_B = @. exp(-im*diagonal_H_B*t/N)

    for _ in 1:N
        psi .*= diagonal_U_A
        psi = fft(psi)
        psi .*= diagonal_U_B
        psi = ifft(psi)
    end
return psi
end

your case

However for KrylovKit you just require the action of H itself and you can actually compute this in a similar way but exact:

H |\psi\rangle = H_A |\psi\rangle + H_B|\psi\rangle

You compute each part of the sum the same way that you would compute it in the split-operator method. So say |\psi\rangle is given in real space where H_A is diagonal and H_B is diagonal in momentum space, then you’d do roughly:

function H(psi)
    diagonal_H_A = ...
    diagonal_H_B = ...

    psi_A = psi .* diagonal_H_A
    psi_momentum_space = fft(psi)
    psi_momentum_space .*= diagonal_H_B
    psi_B = ifft(psi_momentum_space)
    return psi_A + psi_B
end

Then simply use KrylovKit.eigsolve(H, psi_0) done. All you need to do really is construct the diagonals of H_A and H_B. I am not sure whether there are packages that can help you out-of-the-box. Perhaps WignerSymbols.jl is useful if you need Clebsch-Gordon coefficients.

9

Thanks for the insight. The only tricky bit would be diagonalizing the electron-electron repulsion (1/|\vec{r_1} - \vec{r_2}|) but in principle there are routines to do that so long as I can construct the matrix. I’ll do some more looking into the matrix-construction end.