/*

* Copyright (C) 2013-2014 Ulf Lorenz

*

* This program is free software; you can redistribute it and/or modify

* it under the terms of the GNU General Public License as published by

* the Free Software Foundation; either version 2 of the License, or

* (at your option) any later version.

*

* This program is distributed in the hope that it will be useful,

* but WITHOUT ANY WARRANTY; without even the implied warranty of

* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

* GNU Library General Public License for more details.

*

* You should have received a copy of the GNU General Public License

* along with this program; if not, see <http://www.gnu.org/licenses/>.

*/

#include <gtest/gtest.h>

#include "wavepacket.hpp"

using namespace tensor;

using namespace wavepacket;

/** Acceptance test: Propagate a one-dimensional free particle.

*

* We use a Gaussian function as initial state, which is up to normalization

*

* psi(t=0) = exp( - alpha0 * x^2 ),

*

* and just let it evolve in time. The solution for this problem is known analytically,

* and gives (up to normalization and using atomic units)

*

* psi(t) = exp( - alpha(t) * x^2 + i gamma(t) )

*

* where i is the imaginary number, and

*

* alpha(t) = alpha0 / (1 + 2*i*alpha0*t/m)

* gamma(t) = i/2 * ln(1 + 2*i*alpha0*t/m)

*

* At each time step, we compare the overlap between the propagated wave function

* and the analytic reference.

*

* It turns out that this works pretty well even with relatively imprecise

* propagation. However, as soon as the wave function reaches the grid boundaries,

* the analytical and numerical solutions quickly diverge. That is why we use a

* relatively short propagation time and require high fidelity.

*/

class AcceptanceTest : public ::testing::Test {

public:

/** Helper function: Construct the analytic solution without normalization. */

CTensor makeAnalyticSolution(const RTensor& xgrid, double t)

{

cdouble im{0,1};

// the analytic solution up to normalization

cdouble alphat = this->alpha0 / (1.0 + 2.0*im*this->alpha0*(t)/this->mass);

cdouble gammat = im/2.0 * log(1.0 + 2.0*im*this->alpha0*(t)/this->mass);

return exp(-alphat * xgrid*xgrid + im * gammat);

}

// mass of the particle

double mass = 1.0;

// initial width of the Gaussian

double alpha0 = 1;

// Grid parameters: minimum, maximum, #points

double xmin = -20.0;

double xmax = 20.0;

unsigned int npoints = 256;

// time parameters

double t0 = 0;

double tend = 3.9;

double dt = 1;

};

TEST_F(AcceptanceTest, FreeParticle) {

// Set up the grid and the Hamiltonian

auto grid = std::make_shared<FFTGrid>(this->xmin, this->xmax, this->npoints);

auto representation = Representation::create({grid});

auto hamiltonian = std::make_shared<SchroedingerOperator>(

std::make_shared<CartesianKineticEnergy>(representation, 0, this->mass),

true);

// Set up the initial wave function.

auto builder = std::make_shared<StateBuilder>(representation);

CTensor psi0 = builder->buildProductWavefunction({

std::make_shared<TensorState>(makeAnalyticSolution(grid->getDvrGrid(), 0.0))

});

ASSERT_LT(abs(norm2(psi0) - 1), 1e-10);

// Propagate the previously defined Gaussian with a Runge-Kutta solver.

CTensor psi = psi0;

auto propagator = std::make_shared<OdePropagator_RK5>(hamiltonian, 1e-2, 1e-10);

for (double t = this->t0; t < this->tend; t += this->dt) {

// Propagate to the next time step.

psi = propagator->propagate(psi, t, t+this->dt);

ASSERT_NEAR(1.0, norm2(psi), 1e-6)

<< "Norm was not conserved. " << std::endl;

// check that the overlap with the analytic solution is small enough

// to blame it on numerical errors.

CTensor reference = builder->buildProductWavefunction({

std::make_shared<TensorState>(makeAnalyticSolution(grid->getDvrGrid(), t+this->dt))

});

ASSERT_LE(abs(scprod(reference,psi) - 1.0), 1e-5)

<< "Wave function is not the correct Gaussian; overlap too poor." << std::endl;

}

}