Computational implementation

The results in equations 5.19, 5.20, 5.22 and 5.23 have been written in a form which shows that in general each term can be represented by a real numerical prefactor, integers which are the powers of $\{R, q_{n \ell}, q_{n' \ell'} \}$ and one further integer to signify the presence of one of the terms $\{ {\sin q_{n \ell} R},{\sin q_{n' \ell'} R}, {\cos q_{n \ell} R}, {\cos q_{n' \ell'} R} \}$. When these terms are combined and differentiated by the ${\hat D}_{\ell m}$, the general term also needs integers to represent powers of $\{x_{\alpha \beta}, y_{\alpha \beta}, z_{\alpha \beta}, R_{\alpha \beta}, r_{\alpha}, r_{\beta} \}$. Therefore a general term in the expressions for ${\cal S}_{\alpha, n \ell m ; \beta, n' \ell' m'}$ and ${\cal T}_{\alpha, n \ell m ; \beta, n' \ell' m'}$ could be represented by a data structure consisting of one real variable $g$ and ten integer variables $I_{1-10}$ as follows:

$$g \frac{x_{\alpha \beta}^{I_1} y_{\alpha \beta}^{I_2} z_{\al... ...\rightarrow \{ g,I_1,I_2,I_3,I_4,I_5,I_6,I_7,I_8,I_9,I_{10} \}$$ (5.47)

with the following correspondence between $f(R)$ and $I_{10}$:

$$\qquad f(R) = \left\{1 , {\sin q_{n \ell} R} , {\sin q_{n' \ell'} R} , {\cos q_{n \ell} R} , {\cos q_{n' \ell'} R} \right\}$$

$$* \rightarrow I_{10} = \left\{ 0 , 1 , 2 , 3 , 4 \right\} .$$ (5.48)

A recursive function can be written to manipulate these encoded terms and perform the differentiation by the ${\hat D}_{\ell m}$, which can themselves be generated using the recursion rules for the associated Legendre polynomials. Thus it is straightforward to write a code which starts from equation 5.14 and generates the results up to arbitrary values of $\ell$ for $S_{\alpha \beta}$ and $T_{\alpha \beta}$ for the cases when ${\bf R}_{\alpha \beta} \not= 0$. The results for ${\bf R}_{\alpha \beta} = 0$ are simple enough to be coded within the program which uses this basis. For a given ionic configuration the matrix elements between the basis states can be calculated initially and stored on disk for use during the calculation.

The cost of calculating the analytic matrix elements increases dramatically as higher angular momentum components are included. In general, a much smaller value of ${\ell_{\mathrm{max}}}$ is used than is ``recommended'' by the kinetic energy cut-off. However, these basis functions are being used to describe functions localised in overlapping regions, and in this instance, a degree of ``under-completeness'' is desirable. If the basis functions formed a complete set (up to a given kinetic energy cut-off) in each support region, then a variation which is confined to the overlapping region can be equally described by variations in either region. Symmetric and antisymmetric combinations of these variations can be formed, the antisymmetric variation vanishing and thus leaving the density-matrix invariant. Therefore this superposition results in directions in the parameter space with very small curvature which degrade the efficiency of minimisation algorithms (see section 6.2.3). When working with overlapping support functions, it is therefore better to treat ${\ell_{\mathrm{max}}}$ as a convergence parameter along with, rather than derived from, $E_{\mathrm{cut}}$.