Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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
    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 fftw++-1.13/ -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;

  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
  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
      /// 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
      /// 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");
    //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";
share|improve this question

migrated from Sep 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. 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.


Ok, so it looks like your problem is with the meaning of the index on the Fourier-transformed vector. This is impossible to resolve in a general manner as it will depend on the specific conventions used by your FFT algorithm but the general advice is independent of implementation.

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} $$

share|improve this answer
thank you this helps a lot. I'll do a bit more reading. i need to figure out what fftw does. it doesn't take parameters and just gives out a momentum vector of the same size as the position one. – kirill_igum Oct 8 '13 at 19:48

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.