# Split operator FFT quantum dynamics for a harmonic oscillator

I would like to do a numerical quantum dynamics of a displaced gaussian in harmonic oscillator using split-operator method (see bottom of these notes by Hal Evans for the algorithm).

I have a problem with the step (3). I'm not sure how to represent the momentum term $p^2$. Is it just an index of the loop?

creation of the initial state:

  unsigned int n=1400; ///> rank of the state vector psi(x). (x_0, x_1, ..., x_1399)^T
const double dx=0.01;///> step of the state vector psi(x). dx = x_1-x_0;

array1<Complex> psix(n,align); ///> array1 is the array class from fftw package that is wrapped into fftwpp. psi(x) is the state vector in the position space

double x0 = 4; ///> the mean (center) of the gaussian of the initial state.

for(unsigned int i=0; i < n; i++) { ///> go over all points in the position space
double x = (i-0.5*n)*dx; ///> for each point, calculate x
x+=x0; ///> shift the gaussian back by x0
psix[i]=exp(-pow(x,2)); ///> create the gaussian for the initial state
}


here is the part of the code that is doing the mentioned algorithm.

  /// step 1
for (size_t i = 0; i < n; ++i) {
double x = (i-0.5*n)*dx; ///> get value of x for each element of position space vector
ax[i]=exp(-ii*u(x)*dt*0.5*invhb)*psix[i];///> propagate half of the step using potential energy operator U:  a(x)=e^(-i*u*dt/(2*hb))*s(x,t)
}
/// step 2
Forward.fft(ax,fs);///> do a fft into momentum space.
/// step 3
for (size_t i = 0; i < n; ++i) {
double p=i;///> THE PROBLEM IS HERE!!! , assign the value of the momentum to be used in kinetic energy operator. I just used the omega from http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Definition
bp[i]=exp(-ii*p*p*dt*0.5*invm)*fs[i];///> B(p) ≡ e^( −i p2 Δt / 2m) A(p) ///> act with a kinetic energy operator onto the momentum-space state vector.
}
/// step 4
Backward.fftNormalized(bp,gfs);///> inverse fourier transform back into the position space
/// step 5
for (size_t i = 0; i < n; ++i) {
double x = (i-0.5*n)*dx;///> convert indices of the position-space state vector into values of position
psix[i]=exp(-ii*u(x)*dt*0.5*invhb)*gfs[i];///> propagate with potential energy operator through the other half-step: psi(x, t+Δt)  =  e^(−i U Δt / 2ℏ) B(x)
}


Here is my full c++ code (based on the tutorial to fftwpp wrapper example to fftw):

#include <limits>
#include <vector>
#include <complex>
//#include "gnuplot_i.hpp"
#include "Array.h"
#include "fftw++.h"

// instal fftw and
// Compile with:
// g++ -std=c++11 -I fftw++-1.13 -fopenmp main.cc fftw++-1.13/fftw++.cc -lfftw3 -lfftw3_omp; ./a.out

using namespace std;
using namespace Array;
using namespace fftwpp;

double u(double x) {
return x*x;
}
int main()
{
const Complex ii(0.0,1.0);
const double invhb = 1.0;
const double invm = 1.0;

unsigned int n=1400;
const double dx=0.01;
const double dt=0.000001;
unsigned int np=n/2+1;
size_t align=sizeof(Complex);

array1<Complex> psix(n,align);
array1<Complex> fs(n,align);
array1<Complex> gfs(n,align);

fft1d Forward(n,-1,psix,fs);
fft1d Backward(n,1,fs,gfs);

double x0 = 4;

for(unsigned int i=0; i < n; i++) {
double x = (i-0.5*n)*dx;
x+=x0;
psix[i]=exp(-pow(x,2));
}

ofstream ofs("psix.dat");
for (size_t i = 0; i < n; ++i) {
double x = (i-0.5*n)*dx;
ofs << x << "  " << real(psix[i]) << "  " << imag(psix[i]) << "\n";
}

int tn=10;
//// following http://hep.physics.indiana.edu/~hgevans/p411-p610/material/08_pde_iv/schrodinger.html
array1<Complex> ax(n,align); ///> a(x)=e^(-i*u*dt/(2*hb))*s(x,t)
array1<Complex> bp(n,align); ///> B(p) ≡ e^( −i p2 Δt / 2m) A(p)
for (size_t l = 0; l < 100; ++l) {
for (size_t j = 0; j < tn; ++j) {
/// step 1
for (size_t i = 0; i < n; ++i) {
double x = (i-0.5*n)*dx;
ax[i]=exp(-ii*u(x)*dt*0.5*invhb)*psix[i];///> a(x)=e^(-i*u*dt/(2*hb))*s(x,t)
}
/// step 2
Forward.fft(ax,fs);
/// step 3
for (size_t i = 0; i < n; ++i) {
double p=i;
//double p=i*dx;
bp[i]=exp(-ii*p*p*dt*0.5*invm)*fs[i];///> B(p) ≡ e^( −i p2 Δt / 2m) A(p)
}
/// step 4
Backward.fftNormalized(bp,gfs);
/// step 5
for (size_t i = 0; i < n; ++i) {
double x = (i-0.5*n)*dx;
psix[i]=exp(-ii*u(x)*dt*0.5*invhb)*gfs[i];///>psi(x, t+Δt)  =  e^(−i U Δt / 2ℏ) B(x)
}
}
//// to plot gnuplots on the fly using gnuplot++, uncomment the lines below.
//vector<double> xv,psir;
//for (size_t i = 0; i < n; ++i) {
//  xv.push_back((i-0.5*n)*dx);
//  psir.push_back(real(psix[i]));
//}
//Gnuplot g("lines");
//g.plot_xy(xv,psir);
//std::cin.ignore( std::numeric_limits <std::streamsize> ::max(), '\n' );
}

ofstream offs("fs.dat");
for (size_t i = 0; i < n; ++i) {
double p=i*dx;
offs << p << "  " << real(fs[i]) << "  " << imag(fs[i]) << "\n";
}

ofstream ofgfs("psixt.dat");
for (size_t i = 0; i < n; ++i) {
double x = (i-0.5*n)*dx;
ofgfs << x << "  " << real(psix[i]) << "  " << imag(psix[i]) << "\n";
}
}

-

## migrated from physics.stackexchange.comSep 25 '13 at 17:45

This question came from our site for active researchers, academics and students of physics.

One word about your code, Please add comments (one detailed comment for 5 lines is a good scale). A good code must be quickly read by anybody who has never seen this code. Unfortunately, most developpers (even good ones) make the same error... Your step 3 doesn't seem wrong, however I wonder how works the correspondence between $x=0$, or $p=0$ and arrays with indices from $0$ to $n-1$. There should have a additional method to identify $x=0$ (or $p=0$) with some indice, for instance $(n/2)$ – Trimok Sep 25 '13 at 17:03
@Trimok that's what i was thinking about the indices but I think it's one-to-one correspondence. en.wikipedia.org/wiki/Discrete_Fourier_transform#Definition shows that there are $n$ of $x_i$ and $n$ of their fourier transforms $X_i$. – kirill_igum Sep 25 '13 at 17:53

I have not looked at your code and this is not the place for it (though Computational Science might be), but your question looks simple in essence.

You have already transformed your wavefunction $\psi(x)$ into a momentum space representation $A(p)$, and you want to propagate it using the kinetic energy part of the hamiltonian, $p^2/2m$. You are then asked to

Evaluate $$B(p) ≡ e^{−i p^2 Δt / 2m} A(p)$$ in momentum space.

Here $A(p)$ and $B(p)$ are Hilbert space vectors, in a momentum representation, and the exponential is a matrix which is also diagonal in that representation.

In any numerical implementation, of course, $A(p)$ will be discretized to a vector with a finite number of complex entries representing a set of amplitudes $(A(p_i))_{i=1}^N$. You then need to multiply each of these amplitudes by the relevant phase, to leave $$(B(p_i))_{i=1}^N=(e^{−i p_i^2 Δt / 2m}A(p_i))_{i=1}^N.$$

I hope this makes it clear what you need to do. If you have a specific question of how to implement this you should ask at Computational Science, or provide a detailed outline of what your code is doing and which part of the code implements the troublesome step.

In most implementations (and yours looks like one of these) you will change the continuous function $A(x)$ for a vector a(i) which has a discrete index i=0,...,n and which samples the continuous function in a finite range and with a finite resolution, a(i)$=A(-x_0+i\delta x)$. You are then doing a discrete Fourier transform which is of the form $$b(j)=\sum_{i=0}^n e^{2\pi \text i \,ij/n^2}a(i),$$ but unfortunately this is relatively different to your original intention of doing a continuous transform as $$B(p)=\int_{-\infty}^\infty\text dx e^{\text i p x/\hbar} A(x).$$ As a result of this, the computed DFT b(j) will in general be a reasonable approximation of the continuous Fourier transform, but the relationship between them can vary depending on conventions. The momentum resolution and the momentum range you'll have available will, in general, be inversely proportional to the position range and position resolution at which you sampled your original function.
There is one further provision which is very important. The DFT is periodic in its transformed variable, and this means that most DFT implementations will have the negative frequencies in the second half of the range. This is what follows from the definition above, though of course you need to check with your implementation. As a result, one typically has to perform some sort of "folding" of the domain, and the interpretation of b(i) will typically be of the form $$b(j)=\begin{cases} \vphantom{\sum}B(j\delta p), &0\leq j\leq n/2 \\ \vphantom{\sum}B((j-n)\delta p),& n/2\leq j \leq n.\end{cases}$$