# Computer simulation of Schrödinger equation [duplicate]

I am looking for a computer program which simulates the Schrödinger equation (say for a single particle) in two dimensions and for potentials and initial states specified by the user. Typical situations I have in mind are the double slit experiment, harmonic oscillator, Aharonov–Bohm, tunneling effect etc. Unfortunately I was only able to find the one dimensional case on the web.

So are there some 2d programs available?

PS: I don't think my question is a duplicate of this one. I only ask about a single particle, not about thousands. On the other hand I want freedom in my choice of potential. Also I want a numeric solution of a PDE, while the other question is about (as far as I understand, which is not ver far) a Monte Carlo simulation. Anyway these are completely different requirements.

-

## marked as duplicate by Kyle Kanos, John Rennie, Qmechanic♦Apr 2 '14 at 17:27

Are you looking for a program that will do all the work for you, and just show results, or are you looking for a programming environment to write it yourself in? –  Nuclear_Wizard Apr 2 '14 at 10:31
Possibly related question: physics.stackexchange.com/q/18969 –  mpv Apr 2 '14 at 10:31
Related: Are there simple ways to numerically solve the time-dependent Schödinger equation?, now at Computational Science. –  Emilio Pisanty Apr 2 '14 at 10:51
@nuke wiz I am looking for both. A program which does all the work for me and which is easy to modify. –  Jan Weidner Apr 2 '14 at 11:51
This question appears to be off-topic because it is asking for a program (c.f. this Meta post about on/off-topic computational questions) and not about physics. –  Kyle Kanos Apr 2 '14 at 13:54

This probably isn't exactly what you're looking for, but if you're looking for the time-independent bound states of a system, the Fourier grid Hamiltonian method may be applicable. Here is an application of it to the following strange-looking potential well:

Here are a few low-energy bound states:

And here are some of the high-energy states:

I used a fairly symmetric-looking potential surface, but you can do arbitrary shapes or squiggly lines if you really wanted to.

Here is Mathematica code used to draw the potential and compute the bound states. Beware that in 2D, just using 81 grid points in each dimension means that you have to diagonalize a $6561\times 6561$ matrix, which can take several minutes. The following defines the spatial axes and constructs the kinetic energy matrix:

\$HistoryLength = 0;
hue = Compile[{{z, _Complex}}, {(1.0 Arg[-z] + \[Pi])/(2 \[Pi]),
Exp[1 - Max[Abs[z], 1]], Min[Abs[z], 1]},
CompilationTarget -> "C", RuntimeAttributes -> {Listable}];
vec[A_] := Flatten[Transpose[A, Reverse[Range[ArrayDepth[A]]]]];
nx1 = 8;
np1 = 8;
nn1 = nx1 np1;
L1 = 8;
\[CapitalDelta]x1 = L1/nx1;
\[CapitalDelta]p1 = 2 \[Pi]/\[CapitalDelta]x1;
\[Alpha]1 = \[CapitalDelta]p1/(2 \[CapitalDelta]x1);
pMax1 = \[CapitalDelta]p1 np1/2;
axis1 = (Range[nn1] - 1/2) L1/nn1;
TRow1 = If[EvenQ[nn1],
1/2 Append[
Table[(2.0 pMax1^2)/nn1^2 (-1)^j/Sin[\[Pi] j/nn1]^2, {j,
nn1 - 1}], pMax1^2/3 (1 + 2.0/nn1^2)],
1/2 Append[
Table[(2.0 pMax1^2)/nn1^2 ((-1)^j Cos[\[Pi] j/nn1])/
Sin[\[Pi] j/nn1]^2, {j, nn1 - 1}],
pMax1^2/3 (1 + 1.0/nn1^2)]];
T1 = MapIndexed[RotateRight, ConstantArray[TRow1, nn1]];
nx2 = 8;
np2 = 8;
nn2 = nx2 np2;
L2 = 8;
\[CapitalDelta]x2 = L2/nx2;
\[CapitalDelta]p2 = 2 \[Pi]/\[CapitalDelta]x2;
\[Alpha]2 = \[CapitalDelta]p2/(2 \[CapitalDelta]x2);
pMax2 = \[CapitalDelta]p2 np2/2;
axis2 = (Range[nn2] - 1/2) L2/nn2;
TRow2 = If[EvenQ[nn2],
1/2 Append[
Table[(2.0 pMax2^2)/nn2^2 (-1)^j/Sin[\[Pi] j/nn2]^2, {j,
nn2 - 1}], pMax2^2/3 (1 + 2.0/nn2^2)],
1/2 Append[
Table[(2.0 pMax2^2)/nn2^2 ((-1)^j Cos[\[Pi] j/nn2])/
Sin[\[Pi] j/nn2]^2, {j, nn2 - 1}],
pMax2^2/3 (1 + 1.0/nn2^2)]];
T2 = MapIndexed[RotateRight, ConstantArray[TRow2, nn2]];


The following allows you to enter a potential surface and plot it (the code below takes Potential to be the strange crosshair-shaped surface I showed at the beginning, but it can be changed to other shapes):

\[Delta] = 0.3;
Potential[x1_,
x2_] := -600 \[Delta]/(Sqrt[(x1 - L1/2)^2] + \[Delta]) \[Delta]/(
Sqrt[(x2 - L2/2)^2] + \[Delta]) -
300 \[Delta]/(
Abs[Sqrt[(x1 - L1/2)^2 + (x2 - L2/2)^2] - 2] + \[Delta]);
Hamiltonian[x1_, p1_, x2_, p2_] := (p1^2 + p2^2)/2 + Potential[x1, x2];
DensityPlot[Hamiltonian[x1, 0, x2, 0], {x1, 0, L1}, {x2, 0, L2},
PlotRange -> All, PlotPoints -> 90, ColorFunction -> GrayLevel]


The following constructs the Hamiltonian and diagonalizes it:

V = DiagonalMatrix[
vec[Table[
Potential[axis1[[k1]], axis2[[k2]]], {k1, nn1}, {k2, nn2}]]];
T = KroneckerProduct[T2, IdentityMatrix[nn1]] +
KroneckerProduct[IdentityMatrix[nn2], T1];
H = T + V;
{\[Lambda], \[CapitalLambda]} = {#1[[Ordering[#1]]], \
#2[[Ordering[#1]]]} & @@ Eigensystem[H];


And the following plots the bound states:

mat[v_] := Partition[v, nn1]\[Transpose];
interp[m_] :=
8 InverseFourier[
RotateLeft[
8 Dimensions[m]], Floor[Dimensions[m]/2]]];
Manipulate[{\[Lambda][[k]],
If[Abs[\[Lambda][[k - 1]] - \[Lambda][[k]]] <= 0.01 ||
Abs[\[Lambda][[k + 1]] - \[Lambda][[k]]] <= 0.01, "Mixed",
"Pure"],
Image[hue[20 interp[mat[\[CapitalLambda][[k]]]]], ColorSpace -> Hue,
Magnification -> 1]}, {k, 1, 500, 1}]


Unfortunately I don't exactly remember what physical units system the code uses, but I do remember that I wrote it partially based on David Tannor's exposition of the FGH method as outlined in his book Quantum Mechanics: A Time-Dependent Perspective.

-
Fantastic plots and thanks for the code +1. I'm sure a few people will have fun with that. Have you thought about putting it on the Wolfram Demonstrations site? Also, I love your user page, but is there any chance we could get to know some tiny part of the background of someone who clearly has considerable grasp of many fields of physics? –  WetSavannaAnimal aka Rod Vance Apr 3 '14 at 13:16
@DumpsterDoofus although your code works (and I have yet to make some sense of it), I get a warning for the first part: "Compile::nogen: A library could not be generated from the compiled function. >>". –  Ruslan Apr 6 '14 at 4:49

There are some simulation tools available online, but whether they are useful to you depends on the details of your requirements.

Check out the list of quantum simulators here.
Or this one.

-
Thanks! Actually I am looking for something more flexible, where I can play around with different potentials etc. –  Jan Weidner Apr 2 '14 at 11:57