In this appendix we explain in more detail the method for the evaluation of and the method we use to delete points from the grid.
In our program the coefficient vector is ordered with the index running fastest. Thus, the application of [whose elements are defined just below Eq. (7)] is straightforward. is not a function of and j, therefore we store this matrix only once and apply it to the coefficient vector for each combination. Applying to the coefficient vector is not straightforward, because the coefficient vector is ordered with running fastest and not with running fastest. Hence, prior to applying the coefficient vector is reordered with the index running fastest. Afterwards, the resultant vector is reordered back to the original ordering and added to the resultant vector of the operation. Also here, we store the matrix only once. The same procedure is followed for the application of and the rotational terms from Eq. (7), but now we reorder with j running fastest. Of course, actually resorting the vector each time is costly. Therefore, we only sort the vectors once and store the permutations needed to get from one ordering to another (and back again). This procedure is taken from Ref. . Any general purpose routine can be used in this algorithm. We used the HEAPSORT algorithm from Numerical Recipes. After calculating the kinetic energy terms, the potential terms and the rotational terms, we then calculate the Coriolis terms for substate by first getting the coefficient vector for substate from the neighboring processor (using MPI) and subsequently calculating the Coriolis coupling between substates and . Then we do the same thing for the Coriolis coupling with substate . At this point we cycle back until we have completed one time-step.
The calculation of the kinetic/potential terms of the Hamiltonian are the most expensive steps, computationally. The storage of the reduced potential matrix is most expensive memory-wise and we do not want to recalculate this matrix at each time-step, since its calculation is also time-consuming. With the parameters given in Table I, we need approximately 180 Mb to store the matrix. However, a significant number of grid points lie in regions, where the potential is such that the wave function will be zero. Therefore, we have devised a way to get rid of the points in those areas.
The algorithm is simple. We determine, for each combination, whether it lies the asymptotic region of 3 separate atoms or close to the repulsive wall. If either of the two situations is the case, we delete the point from our grid. The way we determine where a point lies on the PES is illustrated by Fig. 6. In this figure we plot the angular dependence of the PES at certain characteristic combinations. Curve 1 lies in the region of three separate atoms. There the potential is always lower than the asymptote and almost flat. For curve 2, which lies close to the repulsive wall, the potential is high for all angles. Thus, we discard points which have certain characteristics.
Using this algorithm only slightly complicates the algorithm for the calculation of . In fact, it only has an impact on the application of and to the coefficient vector, because deleting a point means deleting a row and column from and . However, this is straightforward to implement. As a result of deleting points, one could end up with a row of DVR-points which is cut in two (see Fig. 7 around R=7 a.u.). This means in this case a kinetic energy matrix , which consists of two diagonal blocks corresponding to the two sections of the row of DVR-points and two rectangular off-diagonal block connecting the two sections. At this point we can choose to ignore those off-diagonal blocks or use them in the calculation. We examined both options and conclude that there is not much difference between the two. In the end, we choose to include the off-diagonal blocks in the calculation with the rationale that although we assume that the wave function is zero in the gap between the points, it might still be that there is overlap between the wave function on one side and the second derivative of the wave function on the other side of the gap.