Moving Grids¶

For problems that optimize wavepacket dynamics across large physical distances or where the maximum possible momentum \(k_{\max}\) is expected to be large, it can be advantageous to use moving grids. In this case, instead of having a static grid in both coordinate (\(r\)) and momentum (\(k\)) space (which we will refer to together as phase space), the grid follows the movement of the wavepacket. Not only does this allow one to have dynamics that take the wavepacket outside of the initially defined phase-space boundaries, but it has the additional benefit that the wavefunction has fewer \(r\)-dependent phase components, which makes its representation more efficient, and easier to handle for propagation.

Turning on moving grids¶

To actually make use of moving grids is very straight-forward. First, one must turn on the use of moving grids in the config file by adding the parameter

moveable = true

in the grid section. The alternative spelling movable is also accepted. For example,

grid: n = 1, system = trap
    1 : rmin = -10 rmax = 10, n = 256, moveable = true

This tells the grid that, when doing propagation, it should update the position of the grid at each time step. Once the grid moves, it is likely that operators in coordinate and phase space will have to be refreshed, since we have effectively changed our coordinate system. This is achieved by using a hook function, called upd_hook, which must be passed to any routine that performs propagation, for example prop_step or oct_step.

The `upd_hook` function¶

This function should be implemented as an external procedure in the program. It should follow the interface definition

interface
  subroutine info_hook(psi, grid, ham, ti, work, para)
    use def_mod
    type(psi_t),  intent(in)    :: psi
    type(grid_t), intent(in)    :: grid
    type(ham_t),  intent(in)    :: ham
    integer,      intent(in)    :: ti
    type(work_t), intent(inout) :: work
    type(para_t), intent(in)    :: para
  end subroutine
end interface

This allows a great deal of flexibility in what takes place after each movement of the grid. Typically, we will simply want to refresh our operators on the modified grid by calling init again. In this case, there is a built-in function, init_ops, which does just that. Hence prop_step could be called in the following manner:

prop_step(psi, grid, ham, work, para, ti, bw, pulses, upd_hook=init_ops)

For more detail on the exact implementation, we refer the reader to the example moving_grid in the examples folder.

Implementation¶

The implementation of the moving grid can be divided into two parts: the movement of the grid in coordinate space, and the movement in momentum space. Both these transformations are handled differently.

Movement in r¶

The implementation of the moving grid in \(r\) is very simple. At each time step, we calculate the expectation value of the position operator for the wavefunction \(\psi(t)\);

\[\avg{r} = \bra{\psi(t)} r \ket{\psi(t)}\,.\]

This defines the “middle” of the potential, which we assume should correspond to the center of the \(r\)-grid. We then calculate the difference between \(\avg{r}\) and the actual \(r\)-grid center

\[d_r = \avg{r} - \frac{r_{\max} - r_{\min}}{2}\,,\]

where \(r_{\max}\) and \(r_{\min}\) are the boundaries of the current \(r\)-grid.

Some \(r\)-dependent operators, such as potentials, will be defined externally in a data file. Usually, this data is defined as pairs of points, which describe the value of the potential at a given position. In order to preserve the alignment of the grid to this data, we restrict ourselves to only moving the \(r\)-grid by an integer multiple of the grid spacing \(dr\). Hence we round \(d_r\) to the nearest multiple of \(dr\), which gives us \(\tilde{d}\_r\). Dividing by \(dr\), we return an integer, \(i_r\).

We now have the number of grid-spacings \(i_r\) between the current \(r\)-grid position, and its goal position \(\tilde{d}\_r\). We then remove this many points from one side of the wavefunction, and add grid points to the other side (with a value of 0). This effectively translates the wavefunction. To complete the grid transformation, we add \(\tilde{d}\_r\) to each point of the grid. This transformation is depicted in the following figure.

r-grid movement¶

Movement in k¶

Since the \(k\)-grid has only one operator on it, namely the kinetic operator \(\hbar^2 k^2/2 m\), and this operator is analytical, we aren’t constrained to moving the \(k\)-grid by integer multiples of the \(k\)-grid spacing \(dk\). Hence we follow a slightly different procedure.

Any translation in momentum space can be represented as a phase in coordinate space:

\[\begin{split}\begin{align} \hat{U}_{p}(k^\prime) \tilde{\psi}(k) &= \tilde{\psi}(k + k^\prime) \\ \Rightarrow \quad \hat{U}_{p}(k^\prime) \psi(x) &= e^{i k^\prime x} \psi(x)\,, \end{align}\end{split}\]

where \(\tilde{\psi}(k)\) is the wavefunction in momentum space. At the first propagation step, we calculate the expectation value of the momentum operator

\[\avg{p} = \bra{\psi(t)}\hat{p}\ket{\psi(t)}\]

and store this value as \(k_0\). We then multiply the wavefunction in coordinate space such that

\[\psi^\prime(x,t) = e^{-i k_0 x} \psi(x,t)\]

and then set the wavefunction \(\psi(x,t) = \psi^\prime(x,t)\). We then do the propagation, making sure that we apply the momentum operator on the shifted \(k\)-grid; in other words, the momentum operator becomes \(\hat{T}\frac{\hbar^2}{2m}\left(k - k_0\right)^2\).