Comments (11)
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.
from continuumarrays.jl.
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.
from continuumarrays.jl.
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.
from continuumarrays.jl.
Is this basically radial basis functions (RBFs)? I think that fit's in naturally.
from continuumarrays.jl.
Yes it is.
from continuumarrays.jl.
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.
from continuumarrays.jl.
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:
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)
:
axes(A,1) == OneTo(sum(length.(A.args)))
. This is the standard array choice but obviously doesn't work for infinite axes.axes(A,1) == Inclusion(union(axes(a,1).domain, axes(b,1).domain))
(assuminga
andb
haveInclusion
axes). This is equivalent toPiecewiseSpace
in ApproxFun. It's not clear what to do where the axes intersect.axes(A,1) == axes(a,1) ⊗ OneTo(2)
, (assumingaxes(a,1) === axes(b,2)
). This is equivalent toArraySpace
in ApproxFun. One would access entries asA[SVector(0.1,1)]
from continuumarrays.jl.
PS For completeness, the analogue of SumSpace
is Hcat(a, b)
.
from continuumarrays.jl.
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):
from continuumarrays.jl.
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
from continuumarrays.jl.
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.
from continuumarrays.jl.
Related Issues (20)
- simplify macro does not permit templating? HOT 1
- How to handle Linear operators / functionals? HOT 4
- export diagonal?
- Continuum indices HOT 10
- Continous Linear Algebra HOT 1
- Integral of (basis) functions HOT 1
- Expansion short cut HOT 1
- Split out transform from factorize
- Derivative -> Diff? HOT 2
- Inner product between bases on different grids, basis transforms HOT 1
- Support syntax for kernels HOT 14
- Plot quasivector HOT 5
- [FEATURE]: Issue Template
- [FEATURE]: Pull Request Template
- Support (T/T) \ f for expansions HOT 1
- Product of QuasiDiagonals fails
- Support Plotting HOT 10
- Rename or don't export `grid` HOT 1
- Use bases for operations on AbstractQuasiVector
- Computing bounds on a function over an interval HOT 4
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from continuumarrays.jl.