# Schrödinger equation for a hydrogen atom and lack of memory

I'm trying to solve the Schrödinger equation for a hydrogen atom in the Cartesian coordinate system.

This is my code

h = 1; m = 1; V[r_] := -1/Sqrt[0.000001 + r.r]

\[ScriptCapitalL] = -(h^2/(2 m)) \!$$\*SubsuperscriptBox[\(\[Del]$$, $${x, y, z}$$, $$2$$]$$f[x, y, z]$$\) + V[{x, y, z}] f[x, y, z];

{egv1, eigs1} =
With[{d = 10, n = 3},
NDEigensystem[{\[ScriptCapitalL],
DirichletCondition[f[x, y, z] == 0, True]},
f, {x, -d, d}, {y, -d, d}, {z, -d, d}, n,
Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}},
"Eigensystem" -> {"Arnoldi"}}]];

But my PC is always freezing due to lack of memory. Is there any method to overcome it?

The task has a solution when it is correctly set. Due to the finite size of the region, the eigenvalues do not correspond to the expected values for the hydrogen atom.

h = 1; m = 1; V[r_] := -1/Sqrt[ r.r]

\[ScriptCapitalL] = -(h^2/(2 m)) \!$$\*SubsuperscriptBox[\(\[Del]$$, $${x, y, z}$$, $$2$$]$$f[x, y, z]$$\) + V[{x, y, z}] f[x, y, z];
d = 10; n = 3;
A = ImplicitRegion[x^2 + y^2 + z^2 <= d^2, {x, y, z}];

{vals, funs} =
NDEigensystem[{\[ScriptCapitalL],
DirichletCondition[f[x, y, z] == 0, x^2 + y^2 + z^2 == d^2]},
f, {x, y, z} \[Element] A, n] ;
vals

Out[]= {-0.0122755, -0.0124225, -0.0125422}

Table[
ContourPlot[Evaluate[funs[[i]][x, y, 0]], {x, -d, d}, {y, -d, d},
PlotRange -> All, PlotLabel -> vals[[i]], Contours -> 20,
PlotLegends -> Automatic], {i, Length[vals]}]

The solution in the cubic region (the statement of the author) differs from the solution in the ball in that the eigenvalues become positive, which indicates the influence of boundaries.

h = 1; m = 1; V[r_] := -1/Sqrt[ r.r]

\[ScriptCapitalL] = -(h^2/(2 m)) \!$$\*SubsuperscriptBox[\(\[Del]$$, $${x, y, z}$$, $$2$$]$$f[x, y, z]$$\) + V[{x, y, z}] f[x, y, z];
d = 10; n = 3;
A = ImplicitRegion[-d <= x <= d && -d <= y <= d && -d <= z <= d, {x,
y, z}];

{vals, funs} =
NDEigensystem[{\[ScriptCapitalL],
DirichletCondition[f[x, y, z] == 0, True]},
f, {x, y, z} \[Element] A, n] ;

vals

Out]= {0.00203899, 0.00213474, 0.00233661}

Table[
ContourPlot[Evaluate[funs[[i]][x, y, 0]], {x, -d, d}, {y, -d, d},
PlotRange -> All, PlotLabel -> vals[[i]], Contours -> 20,
PlotLegends -> Automatic], {i, Length[vals]}]

• Hi, could you clarify a bit what you mean with :"...the eigenvalues become positive, which indicates the influence of boundaries." I understand that the region has an effect on the eigenvalues are you suggesting the change in sign can be related to the form of the region? Oct 8, 2018 at 7:06
• Wow! Why is your code running for just 5 sec? I don't see you have spiecified any params for the method at all. Why does it work even with potential singularity at 0 ? Oct 8, 2018 at 11:00
• @JamesFlash, this post explains how NDEigensystem works. Oct 8, 2018 at 15:15
• @AlexTrounev, I do not see how that answers my question. Oct 8, 2018 at 15:17
• @user21Take my code for cube, put V = 0 there, you will get vals={0.0370126, 0.0740349, 0.0740351}. Oct 8, 2018 at 15:22

(More an extended comment than an answer.)

Your spacial resolution is a bit too fine for a common PC.

By using the low-level functionalities of "NDSolveFEM", I was able to assemple the stiffness matrix for "MaxCellMeasure" -> 0.01. It is already 2.8 GB large. A single matrix-vector multiplication (needed for Arnoldi's method) requires about 0.17 seconds on my computer.

You have the options to reduce the resolution, to reduce the interpolation order (from 2 to 1), or to look out for a bigger computer to compute it on.

I also have my doubts that Mathematica's implementation of Arnoldi's method is well-adapted for these large matrices, in particular because we have no way to use preconditioners. Maybe one should employ matrix-free methods (that are not supported by Mathematica at the moment) along with multigrid preconditioners. This should be considerably more efficient, in particular, because we can exploit the tensor-product structure of the mesh grid.

• Thank you for answer Oct 7, 2018 at 18:59