More dimensions

[jagot]

This is a follow-up to #14, but more speculatory.

At the moment, I'm working with atoms with spherical symmetry, which allows me to reduce my differential equations to one-dimensional radial problems, with the angular coordinates treated exactly using angular momentum algebra. In the future, I will probably want/need to go to molecules which possess less/other symmetries. There are many different ways of treating molecules, but the dominating by far is linear combinations of (multi-centred) 3d Gaussians (actually polynomials multiplied by Gaussians). It would be interesting to see if we could fit these type of basis functions within the ContinuumArrays framework. If we get this to work, a Hartree–Fock code could conceivably be used for both atoms and molecules without modifications (some tensor magic would need to happen).

I am by no means well-versed on the topic, and would need to read up quite a bit (although I have a lot of good references). @jarvist could possibly be interested in this as well.

[dlfivefifty]

Yes this needs to be fixed in QuasiArrays.jl. Essentially everywhere that specialises on ::Real is a placeholder for a quasi-array to support general axes types (that is, we want to dispatch on eltype.(axes(A)). Some of these restrictions may be dropped completely, however, we need to be careful with getindex which does different things for vector arguments.

[jagot]

If I understand it correctly, the indexing you're referring to is to evaluate the approximated function at various coordinates (which is certainly necessary). However, I have understood the philosophy of basis sets such as Gaussians to mainly be the transformation of a continuous problem to a matrix problem (with emphasis on this); sort of the opposite what I feel we're trying to accomplish with ContinuumArrays.jl.

I will try to clarify: with molecules, you have no hope to cover all of 3N-dimensional space (N being the amount of electrons) with the same level of accuracy everywhere. Instead, you place a bunch of Gaussians on top of every atom in your molecule and expand your functions using these. There are all kinds of tricks to optimize (reduce) the basis sets for given molecular structures.

Quantum chemists call this "basis set calculations", and you have to derive your equations in terms of these, rather than the actual functions you're approximating, since they do not obey simple things like completeness relations (sum |B_i><B_i| != 1) even when you formally go to infinite amount of basis functions, as do e.g. finite-differences. This is why I'm wondering if we could shoe-horn this into QuasiArrays.jl/ContinuumArrays.jl.

I think my crude explanation is neither fully correct, nor very elucidating, apologies for this.

[dlfivefifty]

Is this basically radial basis functions (RBFs)? I think that fit's in naturally.

[jagot]

Yes it is.

[jagot]

With a little help of @mortenpi, I managed to solve the Dirac equation for hydrogen:

┌─────────┬───────────────────────┬───────────────────────┬─────────────────────────┐
│ Orbital │                     E │                 Exact │                       Δ │
├─────────┼───────────────────────┼───────────────────────┼─────────────────────────┤
│  1s_1/2 │   -0.5000066565965436 │   -0.5000066565953603 │ -1.1833867219479544e-12 │
│  2s_1/2 │  -0.12500208018931686 │  -0.12500208018900594 │  -3.109179580462751e-13 │
│  2p_1/2 │  -0.12500208018919284 │  -0.12500208018900594 │  -1.869060461956451e-13 │
│  2p_3/2 │  -0.12500041602899187 │  -0.12500041602834244 │  -6.494249582544853e-13 │
│  3p_1/2 │  -0.05555629521212829 │  -0.05555629517766647 │  -3.446182228472594e-11 │
│  3s_1/2 │  -0.05555629500066889 │  -0.05555629517766647 │  1.7699758325662174e-10 │
│  3p_3/2 │ -0.055555802242101715 │ -0.055555802093294915 │ -1.4880680021533976e-10 │
│  3d_3/2 │  -0.05555580209449995 │ -0.055555802093294915 │  -1.205036070928145e-12 │
│  4p_1/2 │  -0.03125033815828554 │  -0.03125033803007682 │ -1.2820872141716677e-10 │
│  4s_1/2 │  -0.03125033799437166 │  -0.03125033803007682 │   3.570516105000365e-11 │
│  4p_3/2 │ -0.031250130247878094 │  -0.03125013001044863 │  -2.374294649776232e-10 │
│  4d_3/2 │  -0.03125013001317456 │  -0.03125013001044863 │  -2.725930592362147e-12 │
│  5p_1/2 │ -0.020000181202215195 │ -0.020000181059003808 │ -1.4321138719353144e-10 │
│  5s_1/2 │  -0.02000018113216684 │ -0.020000181059003808 │  -7.316303118898304e-11 │
│  5p_3/2 │ -0.020000074717540672 │ -0.020000074549898272 │ -1.6764239996192032e-10 │
│  5d_3/2 │ -0.020000074512395716 │ -0.020000074549898272 │   3.750255661572055e-11 │
│  6p_1/2 │ -0.013888996756017014 │ -0.013888996749301441 │  -6.715572542503878e-12 │
│  6p_3/2 │ -0.013888935164265204 │ -0.013888935111026512 │ -5.3238691233303825e-11 │
│  6d_3/2 │  -0.01388893511091005 │ -0.013888935111026512 │  1.1646239528317892e-13 │
│  7d_3/2 │ -0.010204112191358305 │ -0.010204112128121778 │  -6.323652712580952e-11 │
└─────────┴───────────────────────┴───────────────────────┴─────────────────────────┘

This is a slightly peculiar equation, since its solutions are vector-valued functions of r, with the vector components customarily referred to as the large and small components, respectively. Since the Dirac Hamiltonian possesses a discrete spectrum (the one in the table above), a positive continuous spectrum, and a negative continuous spectrum, this leads to numeric issues, with spurious eigenstates (due to the operator not being bounded from below, &c). A common way to solve this is to expand the upper and lower components in different basis sets, with e.g. BSplinesQuasi.jl it is possible to use different orders.

What is the correct way to represent the components in an array? For testing, I simply store them as [upper, lower] where upper and lower are column vectors of expansion coefficients, and the Hamiltonian and the basis overlap matrix as BandedBlockBandedMatrices:

This works nicely for B-splines and finite-differences, whose operators are banded, but not for finite-elements, whose operators themselves are BlockBandedMatrices.

Next complication comes when I want multiple functions stored in the same array; for the non-relativistic Schrödinger case, I usually store them as columns: [a b c ...]. I guess the natural extension would be either

• [au bu cu ...; al bl cl ...], or
• a 3d array, M[:,:,1] = [au bu cu ...], M[:,:,2] = [al bl cl ...]

The first idea is easy to implement (and probably efficient), but works only for banded operators, and does not generalize to more dimensions. The second idea is more general, but how do I store operators? I could simply implement a meta-matrix,

struct MetaMatrix{U,UL,LU,L}
   A::U
   B::UL
   C::LU
   D::L
end

and implement mul! &c for that, but that feels a bit ad-hoc.

[dlfivefifty]

In ApproxFun.jl I handled this by interlacing the coefficients by blocks: interlacing rows/columns of banded matrices results in a banded matrix, as does banded-block-banded matrices. I started playing with adding this to LazyArrays.jl as a lazy Interlace array:

JuliaArrays/LazyArrays.jl#41

But in the context of quasi-arrays, some more thought is needed for both concatenation and interlacing. For concatenation, there are multiple options for the axes(A) where A = Vcat(a,b):

1. axes(A,1) == OneTo(sum(length.(A.args))). This is the standard array choice but obviously doesn't work for infinite axes.
2. axes(A,1) == Inclusion(union(axes(a,1).domain, axes(b,1).domain)) (assuming a and b have Inclusion axes). This is equivalent to PiecewiseSpace in ApproxFun. It's not clear what to do where the axes intersect.
3. axes(A,1) == axes(a,1) ⊗ OneTo(2), (assuming axes(a,1) === axes(b,2)). This is equivalent to ArraySpace in ApproxFun. One would access entries as A[SVector(0.1,1)]

[dlfivefifty]

PS For completeness, the analogue of SumSpace is Hcat(a, b).

[jagot]

I see. It seems 1) is basically what I did above, 2) will not work in this case (I think?), and 3) is the most general one, and which is where I want to go.

An operator O acting on the expansion coefficients of such a system will have axes(O) == (axes(a,1) ⊗ OneTo(2),axes(a,1) ⊗ OneTo(2)), I guess? I.e. it is a rank-4 tensor.

[jagot]

I was reading a bit on how to implement the full time-dependent Schrödinger equation for helium, which has two electrons, corresponding to 6-dimensional space. This quickly becomes intractable with growing sizes of the axes, but the nature of the problem is that the electrons are usually never "far away" from the origin simultaneously. It would be cool if we could somehow fit in "Hyberpolic Reduced Tensor-Product Bases", as described by e.g. Lubich (2008):

[dlfivefifty]

This sounds like it’s about tensor Fourier. It wouldn’t be relevant for any basis built from spherical harmonics.

In terms of how it fits in, it’s all about different ways of turning a matrix into a vector. In ApproxFun this was achieved by different tensorizer types, but it was a bit of a hack. We’ll have to address this when support for multiple dimensions is added, hopefully it will be as easy as adding a special tensirizer type

[jagot]

This sounds like it’s about tensor Fourier. It wouldn’t be relevant for any basis built from spherical harmonics.

This is for the radial coordinates of the two electrons, i.e. r_1 and r_2, where the full tensor product basis becomes too large. The angular coordinates, θ_1, θ_2, and φ_12 are treated using bipolar spherical harmonics, they present no problem.

EDIT: I assume you refer to the fact that in e.g. Lubich's book, it's used for tensor products of Hermite functions, which I agree are not suited for spherical bases. My idea was to simply reuse the hyperbolic reduction of the tensor space, but in the radial coordinates. Most publications I've seen go for the full tensor basis, except one, which treats the electrons asymmetrically, i.e. different intervals for 0 <= r_1 <= R_1 and 0 <= r_2 <= R_2, where R_1 << R_2. You can do this because of the symmetry of the problem.