Convergence #1: Equally-spaced grids

I recently got an email from a user that had tried out Wavepacket with some highly accurate data, and found that the results did not match the data. While answering the email, my answer turned out to be a bit long. And I figured out that I may have accumulated quite some experience with respect to converging (wavepacket) calculations over the years. Also, these sort of issues come up every now and again, and there have been similar questions in the past. So I wanted to shape this knowledge into some blog posts about convergence, and how to minimise errors in a wavepacket calculation.

While writing this first blog post, it turned out that properly presenting the maths is rather involved. For this reason, in this post, I only discuss a very special problem: What happens if you have an equally-spaced grid and choose the grid parameters incorrectly? The advantage of equally-spaced grids is that there is a simple and intuitive answer to this question. In later posts I would try to consider other grids, where things become a lot less clear, and also consider time-step problems.

Equally-spaced grids

An equally-spaced grid is determined by three parameters: The center of the grid x_0, the interval length L, and the number of points N.
The latter parameter can also be expressed by the grid spacing . The grid points are equally distributed in the interval , including the two end points.

If you sample any wave function on these grid points, you can do a discrete Fourier transformation to get Fourier components also on an equally-spaced grid in the momentum space. The wave vectors range from -k_max to +k_max with (details depend on whether N is even or odd). For this reason, any wave function that is represented exactly at the grid points (the discrete variable representation DVR) can be equivalently written as a sum of discrete Fourier components (the finite basis representation FBR). For more details, see the wiki page on DVR methods.

In particular, it is common nowadays to use the Fourier method in wavepacket calculations. Here, the momentum operator is applied by transforming into the Fourier space with a discrete Fourier transformation, multiplying each Fourier component with the appropriate momentum and transforming back. Using the Fourier method automatically implies periodic boundary conditions, since the plane waves are generally not limited to the interval L.

Now if we assume for a moment that you manage to choose the parameter x_0 wisely, we are left with the two remaining parameters L and N, which lead to different convergence problems.

Effect of too small interval

If the interval L is chosen too narrow for the dynamics that take place, your wavepacket may travel outside of your grid. When this happens, due to the periodic boundary conditions, the wavepacket may enter the grid from the other side.

For a simple example, let us take a free particle, give it some finite momentum and propagate. The Python code using Wavepacket can be found here, and the resulting video can be seen here.

The initial state is a Gaussian with some positive momentum. As expected, it initially broadens and moves with positive velocity. However, as soon as it reaches the grid boundaries, we can see it entering the grid again thanks to the periodic boundary conditions. Also, towards the end of the video, the Gaussian has broadened so much that the right side (roughly corresponding to the components with the larger momentum) overtakes and interferes with the left side (roughly corresponding to components with smaller momentum).
The resulting oscillations are also typical artefacts if the grid is chosen too small.

In practice, you often encounter another pattern. As an example, let us dissociate a molecule with a Morse potential parametrized for an OH bond. To have a clear wavepacket, the molecule gets a pretty strong kick initially. Again, you can find the code and the resulting video.

What we see here is an initial state with enough energy to dissociate (i.e., a wavepacket travelling to infinity). However, as it reaches the end of the grid, due to the periodic boundary conditions, the wavepacket suddenly feels the influence of the Coulomb barrier at small internuclear distances. As a result, the wavepacket does not reenter the grid from the other side, but is reflected at the boundary of the grid. The barrier here is rather steep, which is a feature not well supported by an equally-spaced grid, but that we can discuss further below.

A problem that you may sometimes encounter is that your grid is always too small. Take the Morse oscillator example, and imagine a laser-driven dissociation of a molecule. The dissociation may take some time, so you may need to propagate for large times. But in these large times, the parts of your wavepacket that dissociated first may reach the boundary of your grid, no matter how large you choose it. What to do then?

Absorbing boundary conditions to the rescue

The solution for this problem are absorbing boundary conditions, typically taking the form of negative imaginary potentials (NIPs). These absorb the wavepacket at the locations where they are nonzero, before the wavepacket reaches the grid boundaries. There are, however, two major drawbacks:

• If the potential is too large, the NIP itself can modify the dynamics.
  Consider a case where the NIP at some grid point absorbs effectively the whole wavepacket in a small timestep. This effect is similar to an infinite potential wall, which would also lead to a density of zero at the grid point, and may thus reflect the wavepacket. A consequence is that you should choose the NIP to be as small and smooth as possible, and you may need to converge this as well.
• With the NIP term, the Hamiltonian is no longer self-adjoint.
  This has serious consequences, because some propagation schemes may no longer work. In particular, the expansion in Chebychev polynomials is only correct for a self-adjoint operator, otherwise you loose all the guarantees (error bounds) that this expansion offers. The workaround here is to move the NIP out of the Hamiltonian, and absorb the wavepacket after every propagation step. This is currently not implemented in the C++ version of Wavepacket, though, because it requires exponentiation of operators.

Let us take the Morse oscillator from the previous section as an example. We might say that an OH distance of 8 atomic units (4 Angstroem) is almost dissociated. Also, from studying the dynamics, it takes at least 100 atomic units (2.5 fs) for the shoulder to travel from 8 a.u. to about the end of the grid. For a constant potential, the absorption of the density is roughly (the factor two comes because the absorption acts on the wavefunction, but the density is the square of the wavefunction). If we want to absorb 99% of the density (A = 0.01), we get a value of V_0 approximately 0.025 H. In practice, we might want to use a smooth turn-on (e.g., harmonic potential) that is zero at r = 8 a.u, and about the intended value at r = 10 a.u., that is, in atomic units, .

Applying this absorbing potential to the Morse oscillator demo gives the results shown in this video (the Python script). For clarity, the y-range has been reduced considerably. We see how the wavepacket is heavily absorbed up to about two percent of the total density. You might want to increase the grid or play around with the NIP if this is not enough, but the principle works as the back of the envelope calculation suggested.

Effect of too few grid points

For a given interval length, the equally-spaced grid also describes an equally-spaced grid of points in the Fourier space. The grid in the Fourier space spans roughly the interval [-k_max, k_max] with . (Details differ for example between even and odd N.)

If this grid is not sufficient to contain the wavepacket, you get artefacts known as aliasing. Roughly, the plane wave and give exactly the same function values at the grid in the real space. But the discrete Fourier transformation always maps this oscillation to the lower wave number k. As a consequence, we have the same effect as the real-space periodic boundary conditions: If the wavepacket accumulates too large Fourier components k > k_max (k < -k_max), it "reenters" the grid at the other side of the Fourier plane

As a simple example, let us consider a linear ramp. A Gaussian wavepacket on the linear ramp experiences some broadening, but also accelerates according to classical equations of motion. If the ramp is sufficiently steep, the wavepacket accumulates sufficient momentum to exceed the range in the Fourier plane.

Let us look at an example, see the code here, and the realspace video here. We can see how the initial wavepacket accelerates towards the negative end of the grid. Shortly before reaching the end of the grid, however, something happens and the wavepacket travels back (upwards on the potential ramp!). The effect is easily understood if we look at the wavepacket in the Fourier space. We see how the wavepacket accumulates momentum until it reaches the end of the grid in Fourier space, and reenters on the side with the positive momenta. After the whole Gaussian has traversed this boundary, we have a Gaussian that has a positive momentum and travels the ramp upwards. This slowly decreases the momentum, and the whole process starts anew.

Hence, equally-spaced grids provide a very simple measure for convergence: eyeballing. Always make sure that your wavepacket stays sufficiently far away from the grid boundaries in Fourier space. Now what is "sufficiently far away"? In principle, the exact definition is a bit complicated and related to off-diagonal operator elements in Fourier space (see the explanation of the DVR approximation). However, for most practical cases, "sufficiently far away" is the typical extension of the Fourier transform of the potential.

Takeaway message

Checking the quality of an equally-spaced grid requires observation of the dynamics both in real space and in Fourier space. This should always be done. In particular, if you find indications of parts of your wavepacket hitting the grid boundaries, you need to increase the interval or number of grid points. These boundary effects are, however, easy to understand and spot. For multi-dimensional problems, you should probably monitor the reduced density.

If there are structural problems, and you cannot make the grid large enough, you can always try to use absorbing boundary conditions. In theory, you could also employ them in the Fourier space, although this seems rather uncommon.

So what is my personal takeaway?
For one thing, writing up such an article is surprisingly hard work. On the other hand, it came back to me that the C++ version of Wavepacket has very poor support for exactly this use-case. The next version (that would be 0.3.3) should definitely feature some usability improvements for plotting and convergence monitoring.