Tutorial.2DHarmonicOscillator

Valid for version 0.1.1.

Here, we start doing a few more complex things, such as setting up multidimensional grids, composing operators, and using the Chebychev propagator.

Program

The following code snipplet is taken from Demos/HarmonicOscillator/Gaussian2D/1/demo.cpp minus the comments:

#include <boost/math/constants/constants.hpp>
#include <wavepacket.hpp>

using namespace wavepacket;
using namespace tensor;
using namespace boost::math::double_constants;

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

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

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

    auto propaPrimitive = std::make_shared<ChebychevPrimitive>(hamiltonian,
            hamiltPrimitive->getSpectralMin(), hamiltPrimitive->getSpectralMax());
    auto propagator = std::make_shared<Propagator>(propaPrimitive, pi/20);
    propagator->addObserver(std::make_shared<LoggingOutput>(hamiltPrimitive), 1);

    propagator->propagate(psi0, 0, 40);

    return 0;
}

In the following, we will only discuss what is new compared to the Free particle tutorial.

Skeleton

We have one new include and using directive.

#include <boost/math/constants/constants.hpp>
using namespace boost::math::double_constants;

This is purely for cosmetic reasons. Among the many things that the boost library provides are various special mathematical functions and constants. Here, we include the definition of pi in the variable boost::math::double_constants::pi. Since this is too long for practical purposes, we import the namespace and can access pi as "pi".

You could have just as well written something like

double pi = 3.141;

Representation setup

This time, we set up a two-dimensional grid,

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

Note that we need to define each one-dimensional grid, then just supply both of them when creating the representation. All of which is a straight-forward extension of the one-dimensional free particle case. With this setup, the first degree of freedom (index 0) always refers to the xgrid, while the second degree of freedom (index 1) refers to ygrid.

Operator setup

This time, we need a more complex operator composed of the kinetic energy for both degrees of freedom plus a harmonic potential in each degree of freedom. You can create such complex operators by simply adding simple operators:

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

The parameters for the operator primitives are again, in this order, the representation to use, the degree of freedom to act on, and the parameters for the operator (mass for kinetic energy, mass and angular frequency for the harmonic potential).

Internally, the summation creates a special sum operator primitive which sums the output of its summands. Operators can also be chained ("multiplied") as "A*B", which produces a new product operator primitive. Finally, it is also possible to add complex operators. I.e., the above code could also be written as

auto primitive1 = std::make_shared<CartesianKineticEnergy>(rep, 0, 1);
auto primitive2 = std::make_shared<CartesianKineticEnergy>(rep, 1, 1);
...
auto hamilt1 = std::make_shared<SchroedingerOperator>(primitive1);
auto hamilt2 = std::make_shared<SchroedingerOperator>(primitive2);
...
auto hamiltonian = hamilt1 + hamilt2 + hamilt3 + hamilt4;

This is rather pointless here, as it only adds additional code and slows down the evaluation a bit. However, Liouvillians, such as commutators or Lindblad Liouvillians, are also implemented as complex operators, and occasionally we want to add them (such as in the case of Lindblad decay, where the total Liouvillian is a commutator plus a Lindblad term).

Initial state

The basic setup is identical to the one-dimensional free particle, only with two functors instead of one, because we have two degrees of freedom.

CTensor psi0 = builder->buildProductWavefunction({
        Gaussian(-5, 0, 1),
        Gaussian(-5, -5, 1) });

As mentioned above, the order of the functors matches the order of the grids. So we have a displaced Gaussian in the x-direction, and a displaced Gaussian with non-zero momentum in the y-direction.

Propagation

Here, we have replaced the simple Runge-Kutta propagator primitive by a much faster propagator based on expansion into Chebychev polynomials,

auto propaPrimitive = std::make_shared<ChebychevPrimitive>(hamiltonian,
        hamiltPrimitive->getSpectralMin(), hamiltPrimitive->getSpectralMax());

The Chebychev propagator takes four parameters: The underlying complex operator that defines the system, the minimum and maximum (real-valued!) eigenvalues of the problem, and an optional accuracy. To get the spectrum, most primitives support helper functions, exceptions being operators that are time-dependent or cannot easily access this information.

Admittedly, supplying this spectral information is pretty annoying. However, guessing reasonable spectral bounds is difficult, especially for operator matrices (which is thoroughly supported by WavePacket), so the design decision was made of not guessing at all for the beginning.