[dd20b8]: test / acceptance / FreeParticle.cpp  Maximize  Restore  History

Download this file

112 lines (98 with data), 4.0 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
/*
* 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;
}
}

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks