Tutorial.FreeParticle

Valid for version 0.1.1

In this tutorial, we discuss the simplest possible wavepacket program: A one-dimensional free particle that we propagate for some time steps while producing some standard log output.

For a typical wavepacket propagation, you need to do the following steps:

1. Set up the representation and the grid
2. Set up the required operators
3. Build the initial state
4. Set up the propagator and propagate

Program

The following code propagates the free particle:

#include <wavepacket.hpp>

using namespace wavepacket;
using namespace tensor;

int main() {
    auto xgrid = std::make_shared<FFTGrid>(-20, 20, 128);
    auto rep = Representation::create( {xgrid} );

    auto hamiltPrimitive = std::make_shared<CartesianKineticEnergy>(rep, 0, 1);
    auto hamiltonian = std::make_shared<SchroedingerOperator>(hamiltPrimitive);

    auto builder = std::make_shared<StateBuilder>(rep);
    CTensor psi0 = builder->buildProductWavefunction({ Gaussian(0, 0, 1) });

    auto propaPrimitive = std::make_shared<OdePrimitive_RK5>(hamiltonian);
    auto propagator = std::make_shared<Propagator>(propaPrimitive, 0.5);
    auto observer = std::make_shared<LoggingOutput>(hamiltPrimitive);

    propagator->addObserver(observer, 1);
    propagator->propagate(psi0, 0, 10);

    return 0;
}

It can be found in the directory Demos/FreeParticle/Gaussian_1D/1/demo.cpp.

The simplest way to compile this program is to take the Makefile in the directory, and possibly adjust the variable WAVEPACKET_TOPDIR to the top level directory of WavePacket. Other required parameters are included automatically.

Now let us go through the code step by step.

Skeleton

The basic skeleton of a wavepacket program looks like this:

#include <wavepacket.hpp>

using namespace wavepacket;
using namespace tensor;

int main() {
    return 0;
}

The include directive makes all WavePacket classes and definitions available, while the two "using" directives allow you to skip the declaration of namespaces (by default, you would have to prepend each WavePacket class by "wavepacket::", e.g., "wavepacket::Representation", because it resides in a namespace of the same name). The main function with zero return value is standard C++.

Representation setup

The grid is set up by these two lines:

    auto xgrid = std::make_shared<FFTGrid>(-20, 20, 128);
    auto rep = Representation::create( {xgrid} );

You set up one independent grid/basis expansion along each degree of freedom. Here, we only have one degree of freedom, and we use an equally-spaced grid (/plane wave expansion) in the first line.

Once you have set up all one-dimensional grids, you combine them into a multidimensional grid, represented by the Representation class. Due to some problems in the C++ standard, the Representation class is an exception to the usual std::make_shared setup; instead, you have to call Representation::create with a list of all one-dimensional grids.

One pecularity should be pointed out: The curly braces here denote a list of items, which are separated by commas. This idiom should always be possible when a list of items is required.

Operator setup

The operators are declared by

    auto hamiltPrimitive = std::make_shared<CartesianKineticEnergy>(rep, 0, 1);
    auto hamiltonian = std::make_shared<SchroedingerOperator>(hamiltPrimitive);

There is an important distinction here between primitive operators and complex operators. Primitive operators know how to act on an arbitrary bra or ket vector. As an example, the CartesianKineticEnergy primitive will perform an elementwise multiplication with the kinetic energy in momentum space for the respective degree of freedom, which on a low level is formally identical for wave functions and density operators.

However, wave functions and density operators have different equations of motion. For wavefunctions, we just apply an operator, while for density operators, we need to commute the Hamiltonian with the density operator, or apply primitive operators for example within a Lindblad Liouvillian. This part of the differential equation is implemented by complex operators. In particular, the SchroedingerOperator here just applies the primitive operator to a wave function.

Initial state

Initial states are built by a StateBuilder object

    auto builder = std::make_shared<StateBuilder>(rep);
    CTensor psi0 = builder->buildProductWavefunction({ Gaussian(0, 0, 1) });

As a sidenote, we want to point out how all objects that need information about the grid get a representation pointer as first argument. This is a common pattern with all objects in WavePacket.

The StateBuilder offers several functions to build initial wave functions and density operators. In the simplest case, we build a multidimensional wave function as a product of wave functions along each degree of freedom. This function requires a list of functors that get the DVR grid for the respective degree of freedom and return the corresponding one-dimensional wave function.

Here, we have only one degree of freedom, and construct a Gaussian function object with zero average position and momentum, and an rms width of 1. Note again that formally we use a list of arguments (hence the outer curly brackets).

Note that we do not use the auto keyword whenever we deal with tensors. The background is that the tensor library uses a copy-on-write mechanism that is, briefly put, incompatible with the auto keyword. Another side note: The brackets in "Gaussian(0,0,1)" can be replaced by curly brackets, which is the new C++-11 syntax. We do not use it in the demos, but it may occur in the code.

Propagation

The work of propagating a state is split up between three types of objects:

    auto propaPrimitive = std::make_shared<OdePrimitive_RK5>(hamiltonian);
    auto propagator = std::make_shared<Propagator>(propaPrimitive, 0.5);
    auto observer = std::make_shared<LoggingOutput>(hamiltPrimitive);

    propagator->addObserver(observer, 1);
    propagator->propagate(psi0, 0, 10);

Similarly to the operators, we also split the propagation between a propagation primitive and a propagator. The propagation primitive gets a complex operator (Hamiltonian or Liouvillian) and encapsulates all details to propagate a given state for a single time step forward or backward in time.

Usually, you do not use the primitive directly. Instead, there is the Propagator class, which, given a primitive and a time step, can propagate a state multiple time steps forward and backward. In this specific example, we set up a standard fifth-order Runge-Kutta solver as primitive, and use an elementary step size of 0.5 atomic units. Internally, the Runge-Kutta solver uses an adaptive step size with a much smaller step; the 0.5 is only the step size that is propagated in one go.

One reason to use a propagator and not the primitive directly is the existence of observers. In the example above, a LoggingOutput observer is registered with the propagator to be called at every ("1") time step of 0.5 atomic units. The propagator will then make sure that this observer is called whenever the propagation has proceeded by one natural time step. The LoggingOutput observer is currently a placeholder that displays a few useful expectation values again roughly corresponding to the output of the Matlab version.

The actual propagation (including the calls to the observers) is finally started with the last line from time 0 to time 10 * 0.5 = 5 with the generated initial wave function as input at time 0.