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.