schroedinger.mws
Solution of the Schroedinger equation by the Numerov-method for any given one-dimensional potential
copyright: Markus Mueller, Hanspeter Huber, Institut fuer Physikalische Chemie, Universitaet Basel
Juli 2000
Introduction and Initialization
The renormalized Numerov method (general annotations and abstract of the theory)
The renormalized Numerov method is widely used to solve the one-dimensional eigenvalue problem for a given potential. The potential can either be a model one (e.g. box potential, harmonic potential) or a potential from a fit to quantum chemically calculated energy values. The central idea is to fix the eigenvalue and solve the differential equation. If the resulting function is not an eigenfunction, the deviation serves as a guide how to change the eigenvalue. Here follows an abstract of the theory.
Schroedinger-equation in a.u.: -
*
+ V(x)
(x) = E
(x)
Transformation u"(x) = f(x)u(x) with u(x) = u =
(x) and f(x) = 2m (V(x) - E)
Taylor-series around x: u(x+h) = u + hu' +
u"+
u"'+
u""+
u""'+
u"""+ ...
I) [u(x+h) + u(x-h)] = 2u +
u" +
u"" +
u"""+ ...
Take the second derivative of I) and multiply with -
:
II) -
[u"(x+h) + u"(x-h)] = -
u" -
u"" -
u""" + ...
I) + II) [u(x+h) + u(x-h)] -
[u"(x+h) + u"(x-h)] = 2u +
u" -
u""" + ...
Replace u" by the differential
equation u"=fu (see above): [u(x+h) + u(x-h)] -
[f(x+h) u(x+h) + f(x-h) u(x-h)] = 2u +
fu - ...
Contract to [1 -
f(x+h)]u(x+h) + [1 -
f(x-h)]u(x-h) = (2 + 10 *
f)u - ...
Higher terms in h are neglected and the following abbreviation is used:
T(x) =
f(x)
Recursion formula: [1 - T(x+h)]u(x+h) + [1 - T(x-h)]u(x-h) = [2 + 10 T(x)]u(x)
The following abbreviations
are used:
F(x) = [1 - T(x)]u(x)
und
U(x) =
![[2+10*T(x)]/[1-T(x)]](/view.aspx?SI=3828/schroedinger49.gif)
F(x+h) = U(x) F(x) - F(x-h) and with the abbreviation:
R(x) =
we get the
2-term recursion formula: R(x) = U(x) -
(x-h)
Similarly with
(x) =
we get an anologous decreasing
recursion
(x) = U(x) - 1 /
(x + h)
We start the integration between
![x[0]](/view.aspx?SI=3828/schroedinger56.gif)
und
![x[1]](/view.aspx?SI=3828/schroedinger57.gif)
with the last recursion formula at x1 (R* =
) and continue until R is about 1 , i.e. the derivative of
is about 0. This point is selected as meeting point
![x[M]](/view.aspx?SI=3828/schroedinger60.gif)
. Then we repeat the procedure with the first recursion formula from
![x[0]](/view.aspx?SI=3828/schroedinger61.gif)
(R =
) until
![x[M]](/view.aspx?SI=3828/schroedinger63.gif)
, and we count how often R<0 , which gives the number of nodes. Meet the two functions, i.e. is DE = R * (
![x[M]](/view.aspx?SI=3828/schroedinger64.gif)
+ h) - R(
![x[M]](/view.aspx?SI=3828/schroedinger65.gif)
) = 0, then the energy E is an eigenvalue and the function
calculated from R resp. R* is an eigenfunction, otherwise the number of nodes and DE give criterias, how E has to be changed.
Further information about the mathematical background of the renormalized Numerov method is provided in the following references:
J. M. Blatt, J. Comput. Phys. 1, 382-396 (1967)
B. R. Johnson, J. Chem. Phys. 67, 4086-4093 (1977).
Starting the program
>
restart;
Digits := 20:
>
max_n_dim := 2000:
This variable should normally not be changed. See also "Advanced utilities" below.
speedup:=10:
This variable should normally not be changed. It means that each eigenvalue is first calculated with an array size, which is speedup-times smaller than max_n_dim to accelerate the calculation. The eigenvalue such found is then used as starting value for the subsequent calculation with the full array size
1.)
Select a potential
In the program section "Selectable potentials" you find a selection of potentials.
2.) Copy the potential procedure name
Press the "plus button" of the subsection for the potential you want to use in the calculation. Copy the procedure name (e.g.
harmonic_potential
:= ....... ).
3.) Insert the potential procedure name
Replace the potential name in the next line on the right of the equal sign (default is: box_potential) by your potential.
>
selected_potential:=box_potential;
4.) Select the orbitals (eigenfucntion) you want to calculate
In the following lines you find the two variables
from_orbital := 1:
and
to_orbital := 5:
Enter your desired orbital range.
>
from_orbital := 1;
to_orbital := 5;
5.) Execute the program
In the "Edit" menu select "Execute" -> "Worksheet".
Advanced utilities
Improvement of the calculation accuracy
The parameter
max_n_dim
("Globals and arrays") determines the array width. Increasing
max_n_dim
(default 2000) leads to more accurate energy levels and higher resolution in the output display. Notice: The greater
max_n_dim
the longer the calculation duration - find a compromise! Our default values are set for teaching, not research! Be aware that an accurate calculation might also need a larger xscale (see below), especially for higher orbitals, because the potentials go to infinity at the end of the array.
Amplitude of the wavefunction
For more informative displays it is sometimes useful to vary the wavefunction amplitude. In every potential procedure you can find the parameter
func_modulator
. Increasing the value of this parameter will increase the amplitude of the wavefunction.
How can I calculate energy levels of my own potential?
In the section
Selectable potentials
a template potential named "
My potential
" is provided at the end. A steady function can easily be implemented by this potential. The constants at the beginning have to be initialised. The function with the variable
x
is inserted at the line
Acquisition[Potential,k] :=
. For non-steady potentials use the potential "Topf (7)" as template. At this point you are already able to start a calculation, but you have to accept that the output is not very informative due to some scale and display parameters that have to be altered.
Local parameters:
Do not touch the parameters k,x and xshift. Just add your own parameters necessary to describe the potential.
Reduced mass:
Calculation of realistic problems (molecule spectroscopy) need finite masses for additional parameters. Usually in vibrational spectroscopy the
reduced mass
is used to describe the system quantitative:
mue := (m1 * m2) /( m1 + m2)
. Use for the mass units atomic units
a.u.
as the energy output of this program is also in
a.u.
->
Microhartree
(for physically meaningful cases like Ne-dimer, H2 etc.)
Display parameters (global):
xscale: default value =1
With the definition "x = (i/max_n_dim -xshift) * xscale" the potential is plotted from x = 0 to x = 1 with the values xshift=0 and xscale=1. For realistic problems one has often to scale the x-axis. For example the Neon2 potential would be plotted and calculated between 0 and 1 a
o
only with xscale=1. For Neon2 xscale is set to 10 and it was necessary to adjust xshift to -0.4 to display all relevant features of the potential in the array.
Attention:
The i-variable of the x-axis is not the distance in Bohr but the array index. Use the formula given at the beginning of this section to calculate the distance.
xshift: default value = 1
Altering the value of xshift will result in shifting the potential in the array boundaries. The harmonic potential for example is shifted by xshift = 0.5. Be careful using this option because shifting important potential features (potential well) out of the array boundaries will result in wrong energy levels. In addition xshift = 0 might produce problems, if the potential goes to infinity at zero.
focus: default value = 1 or aequivalent: round(max_n_dim/max_n_dim)
The parameter focus works as zoom. If you have optimized xshift and xscale and the potential well is displayed too small, you can zoom in by increasing the value of the focus parameter. Altering the focus parameter has no influence on calculation and results because it is only display specific. If the potential well (e.g. Neon2) is not clearly visible, focus according the following recipe:
1.) Check at which array index the left potential flank approximately crosses the x-axis.
2.) Replace the nominator in the default expression"
round(max_n_dim/max_n_dim)"
in such a way that the term becomes a little smaller than the array index determined in the first step.
3.) If the zoom effect is too small, use a smaller number for the nominator.
Notice: It is also possible to insert directly a number for the focus parameter, but then you always have to readjust the focus if you increase accuracy by increasing array resolution (max_n_dim).
Other parameters (global):
func_modulator:
see advanced utilities
potentialscaling:
The parameter
potentialscaling
is a constant factor multiplying the whole potential function. Altering its value will influence the potential energy (y-axis). The energy in the output plot is given in Microhartree by default for physically meaningful cases (Ne-dimer, H2 etc.). If other energy units are recommended this can be done by adjusting the value of the
potentialscaling
parameter.
Globals and arrays
For the input variables see "Introduction and Initialization" -> "Starting the program".
Array dimension of acquisition array and selection of orbitals being calculated.
>
>
numb_orbitals := to_orbital - from_orbital +1:
wavefunction_index := 1:
>
Acquisition := array(0..(2*numb_orbitals),0..max_n_dim):
Stores the orbitals in Acquisition[1..numb_orbitals,i] and lines with the eigenvalues in Acquisition[num_orbitals+1..2*numb_orbitals,i] for the output-graph
Potential := 0:
Acquisition[0,i]=Acquisition[Potential,i] stores the potential function numerically
Epsilon -->Numerov accuracy (error); mue --> default reduced mass
>
epsilon := 1.0e-8:
This accuracy should be good enough for the implemented potentials; for your own potential you might have to adjust it.
>
mue := 1:
Default value for the reduced mass; individual values are set in the potentials-procedures for physically meaningful potentials like Ne2, H2 etc.
Selectable potentials
Box (1)
Particle in a box with infinite high potentials at the walls and zero potential in between
>
box_potential := proc(max_n,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k;
focus := 1:
xscale := 1:
func_modulator := 5:
potentialscaling := 1:
for k from 1 to max_n do
Acquisition[Potential,k] := 0;
od:
>
end:
Harmonic (2)
Particle with a potential that is proportional to its local divergence to normal position.
>
harmonic_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local force_constant,xshift,const,x,k;
force_constant := 4000;
xshift := 0.5;
xscale := 1;
focus := 1;
potentialscaling := 1;
func_modulator := 20;
const := 0.5 * force_constant;
>
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift)*eval(xscale);
Acquisition[Potential,k] := const * x^2;
od:
>
end:
Double I (3)
>
double1_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,b,c,sqr_relative_length,xshift,x;
b := 15000;
c := 2400;
xshift := 0.502;
xscale := 1;
focus := 1;
func_modulator := 20;
potentialscaling := 1;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
sqr_relative_length := x**2;
Acquisition[Potential,k] := b * sqr_relative_length**2 - c * sqr_relative_length +100;
od:
end:
Neon (4)
>
neon_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,Invsqr,Invhoch4,Invhoch6,x,xshift,a1,a2,a3,a4,a5,a6;
a1 := 148.735474;
a2 := 2.3520917;
a3 := 3.749754e4;
a4 := 5.160348e3;
a5 := 8.3724779;
a6 := 6.8028878;
mue := 18395.667;
xscale := 10.0;
xshift := -0.4;
potentialscaling := 1.0e6;
func_modulator := 20;
focus := round(max_n_dim/8);
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
Invsqr := 1 / x**2;
Invhoch4 := Invsqr**2;
Invhoch6 := Invhoch4 * Invsqr;
Acquisition[Potential,k] := potentialscaling * (a1 * exp(-a2 * x) + a3 * (Invhoch6)**2 - a4 * Invhoch4 * Invhoch6 - a5 * Invsqr * Invhoch6 - a6 * Invhoch6);
od:
end:
Morse (5)
>
morse_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,x,D,b,tmpexp,xshift;
D := 5.716;
b := 0.1519399;
xshift := 0.2;
xscale := 40;
focus := round(max_n_dim/18);
potentialscaling := 1;
func_modulator := 1;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
tmpexp := exp(-b * x);
Acquisition[Potential,k] := D * ((tmpexp**2) - 2 * tmpexp);
od:
end:
Invexp (6)
>
invexp_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,x,a,b,xshift;
a := 1000;
b := 5;
xshift := 0.5;
xscale := 1;
focus := round(max_n_dim/max_n_dim);
potentialscaling := 1;
func_modulator := 50;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
Acquisition[Potential,k] := -a * exp(-abs(b * x));
od:
end:
Topf (7)
>
topf_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,x,side,height;
side := 0.2;
height := 200;
focus := round(max_n_dim/max_n_dim);
potentialscaling := 1;
func_modulator := 10;
xscale := 1;
for k from 1 to max_n_dim do
x := (k / max_n_dim)*eval(xscale);
if (x < side) or (x > 1.0 - side) then
Acquisition[Potential,k] := height;
else
Acquisition[Potential,k] := 0;
fi;
od:
end:
Double II (8)
>
double2_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,x,D,w,xo,xshift;
D := 400;
w := 50;
xo := 0.15;
xshift := 0.5;
xscale := 1;
focus := round(max_n_dim/max_n_dim);
func_modulator := 20;
potentialscaling := 1;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
Acquisition[Potential,k] := -D * (exp(-w * (x - xo)**2) + exp(-w * (x + xo)**2));
od:
end:
Wall (9)
>
wall_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,x,width_of_half_wall,height,wo,xshift;
width_of_half_wall := 0.02;
height := 300;
wo := 0.5;
xshift := 0.0;
xscale := 1;
focus := round(max_n_dim/max_n_dim);
func_modulator := 10;
potentialscaling := 1;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
if (x > wo - width_of_half_wall) and (x < wo + width_of_half_wall) then
Acquisition[Potential,k] := height;
else
Acquisition[Potential,k] := 0;
fi;
od:
end:
H2Plus (10)
>
H2Plus_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,x,D,b,c,tmpexp,xshift;
D := 0.109105;
b := 0.707738;
c := -0.379001;
mue:= 919.783;
xshift := 0.09819155;
xscale := 20;
focus := round(max_n_dim/20);
func_modulator := 5000;
potentialscaling := 1.0e6;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
tmpexp := exp(-b * x);
Acquisition[Potential,k] := potentialscaling * (D * ((tmpexp**2) - 2 * tmpexp + c * b * x * ((b * x)**2) * (tmpexp)**2));
od:
end:
HDPlus (11)
>
HDPlus_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,x,D,b,c,tmpexp,xshift;
D := 0.109105;
b := 0.707738;
c := -0.379001;
mue := 1226.37733;
xshift := 0.09819155;
xscale := 20;
focus := round(max_n_dim/20);
func_modulator := 5000;
potentialscaling := 1.0e6;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
tmpexp := exp(-b * x);
Acquisition[Potential,k] := potentialscaling * (D * ((tmpexp**2) - 2 * tmpexp + c * b * x * ((b * x)**2) * (tmpexp)**2));
od:
end:
D2Plus (12)
>
D2Plus_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,x,D,b,c,tmpexp,xshift;
D := 0.109105;
b := 0.707738;
c := -0.379001;
mue := 1839.566;
xshift := 0.09819155;
xscale := 20;
focus := round(max_n_dim/20);
func_modulator := 5000;
potentialscaling := 1.0e6;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
tmpexp := exp(-b * x);
Acquisition[Potential,k] := potentialscaling * (D * ((tmpexp**2) - 2 * tmpexp + c * b * x * ((b * x)**2) * (tmpexp)**2));
od:
end:
Periodic, e.g. metal (13)
>
periodic_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,p,x,D,w,xo,xshift,half_well_numb,Doo,xoo,woo;
D := 400;
w := 50;
xo := 0.15;
half_well_numb := 3;
xshift := 0.5;
xscale := 1;
focus := round(max_n_dim/max_n_dim);
func_modulator := 20;
potentialscaling:= 1;
Doo := 1.4 * D;
xoo := 2 * xo / sqrt(half_well_numb);
woo := 0.08 * w * half_well_numb;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
Acquisition[Potential,k] := 0;
for p from 1 to half_well_numb do
Acquisition[Potential,k] := evalf(Acquisition[Potential,k] - Doo * (exp(-woo * abs(x - (p - 0.5) * xoo)) + exp(-woo * abs(x + (p - 0.5) * xoo))));
od:
od:
end:
Neon2
(14)
>
neon2_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,Invsqr,Invhoch4,Invhoch6,x,x2,xshift,a11,a12,a13,a14,a15,a16,a17,a18,a19;
a11 := 1.174594e6;
a12 := 0.355905;
a13 := 1.325805e5;
a14 := 0.214065;
a15 := 976.69318;
a16 := 0.105389;
a17 := 5.412694e9;
a18 := -3.512808e8;
a19 := -3.964331e6;
mue := 18395.667;
xscale := 10.0;
xshift := -0.4;
focus := round(max_n_dim/8);
func_modulator := 10;
potentialscaling := 1.0e6;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * xscale;
x2 := x**2;
Invsqr := 1 / x2;
Invhoch4 := Invsqr**2;
Invhoch6 := Invhoch4 * Invsqr;
Acquisition[Potential,k] := potentialscaling/1.0e6* (a11 * exp(-a12 * x2) + a13 * exp(-a14 * x2) + a15 * exp(-a16 * x2) + a17 * Invhoch4 * Invhoch6 + a18 * Invsqr * Invhoch6 + a19 * Invhoch6);
od:
end:
NeonHFDB (15)
>
neonHFDB_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,Invsqr,Invhoch4,Invhoch6,x,x2,xshift,a21,a22,a23,a24,a25,a26,a27,a28,a29,Fx;
a21 := 133.8;
a22 := 8.9571795e5;
a23 := 13.86434671;
a24 := 0.12993822;
a25 := 0.24570703;
a26 := 0.53222749;
a27 := 1.21317545;
a28 := 1.36;
a29 := 5.841;
mue := 18395.667;
xscale := 10.0;
xshift := -0.4;
focus := round(max_n_dim/10);
func_modulator := 20;
potentialscaling := 1.0e6;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale) / a29;
x2 := x**2;
Invsqr := 1 / x2;
Invhoch4 := Invsqr**2;
Invhoch6 := Invhoch4 * Invsqr;
if x < a28 then
Fx := exp(-(a28 / x - 1)**2);
else
Fx := 1;
fi;
Acquisition[Potential,k] := potentialscaling/1.0e6*(a21 * (a22 * exp(-a23 * x - a24 * x2) - Fx * (a25 * Invhoch4 * Invhoch6 + a26*Invsqr * Invhoch6 + a27 * Invhoch6)));
od:
end:
NeonEXP6 (16)
>
neonEXP6_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,Invsqr,Invhoch4,Invhoch6,x,x2,xshift,a16,alpha,epsilon,r_asterix;
alpha := 14.66;
epsilon := 115.419;
r_asterix := 5.897836;
mue := 18395.667;
xscale := 10.0;
xshift := -0.4;
focus := round(max_n_dim/9);
func_modulator := 20;
potentialscaling := 1.0e6;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale) / r_asterix;
x2 := x**2;
Invsqr := 1 / x2;
Invhoch4 := Invsqr**2;
Invhoch6 := Invhoch4 * Invsqr;
a16 := alpha - 6;
Acquisition[Potential,k] := potentialscaling/1.0e6*epsilon * (6 / a16 * exp(alpha * (1 - x)) - alpha / a16 * Invhoch6);
od:
end:
ArgonWoon (17)
>
argonWOON_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,x,xshift,wo1,wo2,wo3,wo4,wo5,wo6,expwo54,expwo64;
wo1 := 0.4228;
wo2 := 0.6060;
wo3 := 0.6627;
wo4 := 7.1714;
mue := 36423.46;
xscale := 10;
xshift := -0.4;
focus := round(max_n_dim/4.5);
func_modulator := 50;
potentialscaling := 1.0e6;
wo5 := wo2 * wo3 / ((wo2 + 1.0) * wo1);
wo5 := sqrt(wo5);
wo6 := (wo2 + 1.0) * wo3 / wo2 / wo1;
wo6 := sqrt(wo6);
expwo64 := exp(wo6*wo4);
expwo54 := exp(wo5*wo4);
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
Acquisition[Potential,k] := potentialscaling/1.0e3*(wo1*wo2*exp(wo6*wo4)*exp(-wo6*x)-wo1*(wo2+1.0)*exp(wo5*wo4)*exp(-wo5 * x));
od:
end:
NH3- inversion (18)
>
NH3inv_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
local k,x,xshift,const,a1,a2,a3,a4,a5,a6,a7,a8;
const := -56.3887;
a1 := 0.0000;
a2 := -0.1478;
a3 := 0.0000;
a4 := 0.4248;
a5 := 0.0000;
a6 := -0.2203;
a7 := 0.0000;
a8 := 3.4985e(-2);
Calculation of the reduced mass using: N -> 25,532.71 au and (3*H) -> 5,512.08 au with mue = (m1*m2)/ m1+m2
mue := 4533.4;
xscale := 1.2;
xshift := 0.503;
focus := round(max_n_dim/max_n_dim);
func_modulator := 0.003;
potentialscaling := 1.0;
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * xscale;
Acquisition[Potential,k] := potentialscaling * (const + (a1 * x) + (a2 * x**2) + (a3 * x**3) + (a4 * x**4) + (a5 * x**5) + (a6 * x**6) + (a7 * x**7) + (a8 * x**8));
od;
end:
My Potential (19)
This example potential is the harmonic one and should help the user to implement his own potential by altering the potential function and then adjusting the different parameters step by step until a satisfying output is generated.
>
my_potential := proc(max_n_dim,Acquisition,focus::evaln,mue::evaln,xscale::evaln,potentialscaling::evaln,func_modulator::evaln)
Local parameters:
>
local k,x,xshift,force_constant;
Add here additional parameters of your potential and add their names to the list of local parameters:
force_constant:=4000;
Reduced mass:
>
mue := 1;
Display parameters:
xscale := 1;
xshift := 0.5;
focus := round(max_n_dim/max_n_dim);
Parameter to increase (> 1) or decrease (< 1, > 0) the amplitude of the wavefunction:
func_modulator := 20;
Parameter to alter the energy of the potential by multiplication with a constant -> numerical adjustment to empiric levels:
potentialscaling := 1;
Storage of the potential function into an array (continous -> discrete):
>
for k from 1 to max_n_dim do
x := (k / max_n_dim - xshift) * eval(xscale);
Acquisition[Potential,k] := 0.5*force_constant*x**2;
od;
>
end:
Program procedures
Filling of potential and determination of potential minimum
>
fillpotential := proc(max_n,Acquisition,potential_min::evaln)
>
local infinit,k;
infinit := 100000;
potential_min := infinit:
>
selected_potential(max_n,Acquisition,focus,mue,xscale,potentialscaling,func_modulator);
fills potential into Acquisition[0,i]
>
for k from 1 to max_n do
if Acquisition[Potential,k] < eval(potential_min) then
potential_min := Acquisition[Potential,k];fi;od;
finds potential minimum
>
end:
Startvalues of energy
>
startvalues := proc(max_n,numb_of_wished_nodes,to_orbital,E,Acquisition,E_lower,E_upper,E_lower_keep,point_min)
>
local k, E_upper_start,infinit;
infinit := 100000;
if eval(point_min) then
first trial for each eigenfunction
E_upper_start := -infinit;
>
for k from 1 to max_n do
>
if Acquisition[Potential,k] > E_upper_start then
>
E_upper_start := Acquisition[Potential,k];fi;od;
if Acquisition[Potential,1] = Acquisition[Potential,2] then
E_upper_start := infinit;fi;
if (numb_of_wished_nodes = to_orbital - 1) then
E_lower := potential_min;
E_upper := E_upper_start;
E_lower_keep := eval(E_lower);
>
else E_lower := E_lower_keep;
E_upper := E;fi;
else
if in a first calculation with a smaller array (see speedup above) an eigenvalue E was found that is used as starting value
>
if eval(E) > 0 then
E_lower := 0.9 * E;
E_upper := 1.3 * E;
else
E_lower := 1.1 * E;
E_upper := 0.7 * E;
fi;
>
fi;
end:
Numerov
>
U := proc(n,prefactor,E,Acquisition,T::evaln)
local U_return;
T := eval(prefactor) * (Acquisition[Potential,n] - eval(E));
U_return := (2 + 10 * eval(T)) / (1 - eval(T));
U_return;
end:
>
Numerov := proc(max_n,E,prefactor,nodes::evaln,M,dE::evaln,Acquisition)
>
local k,n,T,nodes_internal;
>
Acquisition[wavefunction_index,max_n] := 1.0e300;
Acquisition[wavefunction_index,1] := 1.0e300;
>
k := max_n;
n := k - 1;
T := 2;
Iteration from the right:
>
while not ((T < 1) and (Acquisition[wavefunction_index,n] <= 1) or (n = 2)) do
n := k - 1;
Acquisition[wavefunction_index,n] := U(n,prefactor,E,Acquisition,T) - 1.0 / Acquisition[wavefunction_index,k];
k := k - 1;
od;
M := n;
nodes_internal := 0;
nodes := nodes_internal;
Iteration from the left:
>
for k from 2 to (eval(M) - 1) do
Acquisition[wavefunction_index,k] := U(k,prefactor,E,Acquisition,T) - 1.0 / Acquisition[wavefunction_index,k - 1];
if (T < 1) and (Acquisition[wavefunction_index,k] < 0) then;
nodes_internal := nodes_internal + 1;
nodes := nodes_internal;
fi;
od;
dE := 1.0 / Acquisition[wavefunction_index,eval(M)+1] - U(M,prefactor,E,Acquisition,T) + 1.0 / Acquisition[wavefunction_index,eval(M)-1]:
>
end:
Sekant- & Bisekant method
>
sekant_method := proc(dE,dEold,E::evaln,Eold::evaln)
local Enew,E_numerator,E_denominator;
E_numerator := eval(E) * dEold - eval(Eold) * dE;
E_denominator := dEold - dE;
if E_denominator = 0 then
Enew := eval(E);
else
>
Enew := E_numerator/E_denominator;
fi;
Eold := eval(E);
E := Enew;
>
end:
>
bisect := proc(nodes,numb_of_wished_nodes,dE,nodeequal::evaln,E::evaln,E_lower::evaln,E_upper::evaln,Eold::evaln)
if (nodes > numb_of_wished_nodes) then
E_upper := eval(E);
elif (nodes < numb_of_wished_nodes) then
E_lower := eval(E);
else
if dE < 0 then
E_lower := eval(E);
else
E_upper := eval(E);
fi;
if (abs(dE) < 1.0) then
nodeequal := true;
fi;
fi;
Eold := eval(E);
E := 0.5 * (eval(E_lower) + eval(E_upper));
end:
FindEigenvalue
>
FindEigenvalue := proc(max_n,numb_of_wished_nodes,to_orbital,point_min,E_lower::evaln,E_lower_keep,E_upper::evaln,prefactor::evaln,E::evaln,epsilon,M::evaln,Acquisition,iteration::evaln)
>
local dE,dEold,Eold,nodeequal,nodes,conv_criterium,Enew,iteration_number;
prefactor := evalf(mue /(6 * potentialscaling * ((max_n - 1) / xscale)**2)):
>
startvalues(max_n,numb_of_wished_nodes,to_orbital,E,Acquisition,E_lower,E_upper,E_lower_keep,point_min);
>
nodeequal := false;
iteration_number := 1;
conv_criterium := 1;
Eold := 20000;
Enew := 0.5 * (eval(E_lower) + eval(E_upper));
>
while epsilon < conv_criterium do
dEold := dE;
conv_criterium := abs((Enew - Eold) / Enew);
Numerov(max_n,Enew,prefactor,nodes,M,dE,Acquisition);
if (nodeequal) and (nodes = numb_of_wished_nodes) then
sekant_method(dE,dEold,Enew,Eold);
else
bisect(nodes,numb_of_wished_nodes,dE,nodeequal,Enew,E_lower,E_upper,Eold);
fi;
iteration_number := iteration_number + 1;
iteration := iteration_number;
od;
E:= Enew:
end:
FindEigenfunction
>
FindEigenfunction := proc(max_n_dim,M,prefactor,E,Acquisition,wavefunction_index)
>
local k,T;
>
Acquisition[wavefunction_index,M] := 1;
>
for k from (M - 1) to 1 by -1 do
Acquisition[wavefunction_index,k] := Acquisition[wavefunction_index,k + 1] / Acquisition[wavefunction_index,k];od;
>
for k from (M + 1) to max_n_dim by 1 do
Acquisition[wavefunction_index,k] := Acquisition[wavefunction_index,k - 1] / Acquisition[wavefunction_index,k];od;
>
for k from 1 to max_n_dim by 1 do
T := prefactor * (Acquisition[Potential,k] - E);
Acquisition[wavefunction_index,k] := Acquisition[wavefunction_index,k] / (1 - T);
od;
>
end:
Acquisition array -> Energy storage & wavefunction modulation
>
Acquisition_modulation := proc(Acquisition,E)
local k;
>
Acquisition[wavefunction_index,0] := E;
for k from 1 to max_n_dim do
Acquisition[wavefunction_index,k] := func_modulator * Acquisition[wavefunction_index,k] + Acquisition[wavefunction_index,0];
Acquisition[(wavefunction_index + numb_orbitals),k] := Acquisition[wavefunction_index,0];
od;
>
end:
Plotting
>
output_plot := proc(Acquisition,numb_orbitals,focus)
local j,color,Acquisition_color;
color[0]:= red:
color[1]:= blue:
color[2]:= green:
color[3]:= orange:
with(plots):
with(plottools):
Acquisition_color := array(0..(2*numb_orbitals));
Acquisition_color[Potential] :=
curve([(seq([k,Acquisition[Potential,k]],k=focus..max_n_dim))],colour=black,style=LINE, thickness=4):
>
for j from 1 to numb_orbitals do
>
Acquisition_color[j] := curve([(seq([k,Acquisition[j,k]],k=1..max_n_dim))],colour=color[j mod 4],style=LINE, thickness=2):
od:
for j from (numb_orbitals + 1) to (2*numb_orbitals) do
Acquisition_color[j] := pointplot({seq([k,Acquisition[j,k]],k=1..max_n_dim)},colour=black,symbol=point):
od:
>
display([seq(Acquisition_color[i],i=0..(2*numb_orbitals))],title="Display of Potential and Orbitals");
you may use as option:
view=[x1..x2,y1..y2] to limit your graph
>
end:
Main program
>
numb_of_wished_nodes := to_orbital - 1:
>
E:=10000:
if from_orbital = 0 or to_orbital = 0 then WARNING("Orbital 0 is not a valid choice !!! - Please interrupt the calculation by pressing the STOP button in the menu and start a new calculation using legitime orbital numbers")
fi;
>
while numb_of_wished_nodes <> from_orbital - 2 do
loop over desired orbitals (eigenfunctions) starting with the energetically highest
max_n := max_n_dim/speedup/10;
point_min := true;
>
while max_n < max_n_dim do
loop over differently sized Acquisition arrays (normally one with size max_n_dim/10 and max_n_dim) to accelerate the calculation (see speedup above)
>
max_n:=max_n*10;
>
fillpotential(max_n,Acquisition,potential_min); FindEigenvalue(max_n,numb_of_wished_nodes,to_orbital,point_min,E_lower,E_lower_keep,E_upper,prefactor,E,epsilon,M,Acquisition,iteration);
point_min := false;
>
od;
end of loop over the differently sized Acquisition arrays
>
FindEigenfunction(max_n_dim,M,prefactor,E,Acquisition,wavefunction_index):
Acquisition_modulation(Acquisition,E);
print(Orbital,numb_of_wished_nodes + 1,"Energy in Microhartree",Acquisition[wavefunction_index,0],"Number of iteration steps",(iteration));
>
numb_of_wished_nodes := numb_of_wished_nodes - 1:
wavefunction_index := wavefunction_index + 1;
>
od:
end of loop over the orbitals
>
output_plot(Acquisition,numb_orbitals,focus);
Warning, the name changecoords has been redefined
![[Maple Plot]](/view.aspx?SI=3828/schroedinger75.gif)
