We choose natural units and set 
, so to recover the ordinary units:
We work in a box of size L and proceed by time steps Δt (in physical units). The length of the box is divided into Nx discrete intervals.
Since the norm of Psi has to be a number, 
implies:
Now for the units of the potential V (in one dimension for brevity):
means
i.e.
hence 
.
We introduce the dimensionful scale 
("electron mass") to make everything dimensionless in the calculations: noting that
,
one translates physical lengths L (and their increments ΔL) to
,
times t to
,
masses m to
,
potentials V to
,
wave numbers k to
and the wave function Ψ to
.
Under the change to dimensionless variables,
becomes then
i.e.
We enforce the requirement that each frequency inserted into the initial wave function be a multiple of a fundamental one: this prevents, in case of periodic boundary conditions, the occurrence of discontinuities.
Physically, we require that eikL=1, which implies
.
The fact that any mode must have a wave number that is a multiple of the above h0 is enforced in the code (by approximating any naively set value).
With the caveat that the value of the potential at its minimum is arbitrary (we customarily set it to zero, however), we calculate the energy as
which in dimensionless units leads to
,
where
We exploit the evolution of the system between successive plotted timesteps: the two wavefunctions (before and after) encode the energy operator, hence (using the fact that the norm of the wavefunction is one)
In this way we do not need to perform additional computation to get the energy, save for a simple scalar product between the two successive wavefunction vectors.
In practice we use the wavefunction before and after a single timestep Δτ, meaning that the results deviates from the true energy in a way that goes to zero as the timestep goes to zero.
In the adimensional units, the energy calculated in this way is given by:
.
Note: this calculation introduces too large a bias in typical configurations, hence it is discarded in favour of the exact one based on the actual definition of energy. The latter, if properly implemented, hitches a ride on the system evolution (at least in the case of the time-dependent potential integrator) and introduces practically no additional overhead.
Periodic and fixed can be chosen; with periodic b.c. one has a little more control over divergencies and spurious high-frequency modes that may arise; however, a proper tuning of the relationship between time and space integration intervals keeps them under control when using the Runge-Kutta scheme described below. In any case the whole-frequency enforcing is performed regardless of the boundary conditions.
The explicit implementation that is used in this code is as follows: calling F the discrete evolution operator, as in
,
the following spatial functions are calculated:
so that for the next update step the wave function has been incremented by
The Runge-Kutta approach above is linear in the wave function, so it can be recast in terms of the action of a matrix H acting on φ to produce the wavefunction for the next timestep.
Moreover, if φ becomes 1+Δφ in each single timestep, then one can perform N updates as one matrix operation after building
,
so that the (N-fold) update reads φ → Uφ. In turn, calling Δτ the "small" timestep, H is built as
,
where
.
Note: In cases where the potential is a function of time, however, all of this pre-processing is not very useful, probably, since most of the computation will need to be performed at each change of the potential (i.e. possibly at each draw operation or even time step).