Internal Diffusion in LJ

We have applied the graph transformation method to study the centre-to-surface atom migration in 55-atom Lennard-Jones cluster. The global potential energy minimum for LJ $_{55}$ is a Mackay icosahedron, which exhibits special stability and `magic number' properties [280,279]. Centre-to-surface and surface-to-centre rates of migration of a tagged atom for this system were considered in previous work [10]. In Reference Wales04 the standard DPS procedure was applied to create and converge an ensemble of paths linking the structure of the global minimum with the tagged atom occupying the central position and structures where tagged atom is placed in sites that lie on fivefold and twofold symmetry axes (Figure 4.15). We have reused this sample in the present work.

**Figure:** Global minimum of the LJ $_{55}$ cluster (shown in stereo). The tagged atom can occupy the central position (green) or one of the two different surface sites (red and blue).
$\begin{psfrags} \centerline{\includegraphics[width=.65\textheight]{markov/arrhenius/LJ55/endpoints/global.ps}} \end{psfrags}$

The sample contained 9907 minima and 19384 transition states. We excluded transition states that facilitate degenerate rearrangements from consideration. For minima interconnected by more than one transition state we added the rate constants in each direction to calculate the branching probabilities. Four digraph representations were created with minimum degrees of 1, 2, 3 and 4 via iterative removal of the nodes with degrees that did not satisfy the requirement. These digraphs will be referred to as digraphs 1, 2, 3 and 4, respectively. The corresponding parameters are summarised in Table 4.1. Since the cost of the GT method does not depend on temperature we also quote CPU timings for the DGT, SGT and SDGT methods for each of these graphs in the last three columns of Table 4.1. Each digraph contained two source nodes labelled and and a single sink. The sink node corresponds to the global minimum with the tagged atom in the centre (Figure 4.15). It is noteworthy that the densities of the graphs corresponding to both our samples (LJ $_{38}$ and LJ $_{55}$ ) are significantly lower than the values predicted for a complete sample [115], which makes the use of sparse-optimised methods even more advantageous. From Table 4.1 it is clear that the SDGT approach is the fastest, as expected; we will use SDGT for all the rate calculations in the rest of this section.

**Table:** Properties of four digraphs corresponding to the LJ $_{55}$ PES sample from an internal diffusion study. **$\vert V\vert$** is the number of nodes; **$\vert E\vert$** is the number of directed edges; **$d_{min}$** , **$\left <d\right >$** and **$d_{max}$** are the minimum, average and maximum degrees, respectively; **$\rho$** is the graph density, defined as a ratio of the number of edges to the maximum possible number of edges; and are the graph radius and diameter, defined as the maximum and minimum node eccentricity, respectively, where the eccentricity of a node is defined as the maximum distance between and any other node. **$\left <l\right >$** is the average distance between nodes. The CPU time, , necessary to transform each graph using the DGT, SGT and SDGT methods is given in seconds for a single 32-bit Intel **$^{\text \textregistered }$** Pentium **$^{\text \textregistered }$** 4 3.00GHz 512Kb cache processor.
$\vert V\vert$	$\vert E\vert$	$d_{min}$	$\left <d\right >$		$d_{max}$	$\rho/10^{-4}$		.		$\left <l\right >$		.		$t_{\rm DGT}$		.	$t_{\rm SGT}$		.	$t_{\rm SDGT}$		.
9843	34871	1	3	.	9	983	3	.	6	10	5	.	71	20	2346	.	1	39	.	6	1	.	36
6603	28392	2	4	.	8	983	6	.	5	9	4	.	86	17	1016	.	1	38	.	9	1	.	33
2192	14172	3	7	.	9	873	29	.	5	4	3	.	63	8	46	.	9	5	.	9	0	.	49
865	7552	4	1	.	9	680	101	.	0	4	3	.	07	7	3	.	1	0	.	8	0	.	12
				.				.				.				.			.			.

For this sample KMC calculations are unfeasible at temperatures lower than about (Here is expressed in the units of $\epsilon/k_B$ ). Already for the average KMC trajectory length is $7.5\times 10^{6}$ (value obtained by averaging over 100 trajectories). In previous work it was therefore necessary to use the DPS formalism, which invokes a steady-state approximation for the intervening minima, to calculate the rate constant at temperatures below [10]. Here we report results that are in direct correspondence with the KMC formulation of the problem, for temperatures as low as .

Figure 4.16 presents Arrhenius plots that we calculated using the SDGT method for this system. The points in the green dataset are the results from seven SDGT calculations at temperatures $T\in\{0.3,0.35,\dots,0.6\}$ conducted for each of the digraphs. The total escape probabilities, $\Sigma_{1}^G$ and $\Sigma_{2}^G$ , calculated for each of the two sources at the end of the calculation deviated from unity by no more than $10^{-5}$ . For higher temperatures and smaller digraphs the deviation was smaller, being on the order of $10^{-10}$ for digraph 4 at , and improving at higher temperatures and/or smaller graph sizes.

At temperatures lower than the probability deviated by more than $10^{-5}$ due to numerical imprecision. This problem was partially caused by the round-off errors in evaluation of terms $1-P_{\alpha,\beta}P_{\beta,\alpha}$ , which increase when $P_{\alpha,\beta}P_{\beta,\alpha}$ approaches unity. These errors can propagate and amplify as the evaluation proceeds. By writing

Terms of the form $1-P_{\alpha,\beta}P_{\beta,\alpha}$ approach zero when nodes $\alpha$ and $\beta$ become `effectively' (i.e. within available precision) disconnected from the rest of the graph. If this condition is encountered in the intermediate stages of the calculation it could also mean that a larger subgraph of the original graph is effectively disconnected. The waiting time for escape if started from a node that belongs to this subgraph tends to infinity. If the probability of getting to such a node from any of the source nodes is close to zero the final answer may still fit into available precision, even though some of the intermediate values cannot. Obtaining the final answer in such cases can be problematic as division-by-zero exceptions may occur.

Another way to alleviate numerical problems at low temperatures is to stop round-off errors from propagation at early stages by renormalising the branching probabilities of affected nodes $\beta \in Adj[\gamma]$ after node $\gamma$ is detached. The corresponding check that the updated probabilities of node $\beta$ add up to unity could be inserted after line 33 of Algorithm B.4 (see Appendix B), and similarly for Algorithm B.3. A version of SDGT method with this modification will be referred to as SDGT2.

Both SDGT1 and SDGT2 have similarly scaling overheads relative to the SDGT method. We did not find any evidence for superiority of one scheme over another. For example, the SDGT calculation performed for digraph 4 at yielded $\mT ^G=\mT ^G_1 W_1 + \mT ^G_2 W_2=6.4\times10^{-18}$ , and precision was lost as both $\Sigma_{1}^G$ and $\Sigma_{2}^G$ were less than $10^{-5}$ . The SDGT1 calculation resulted in $\mT ^G=8.7\times10^{-22}$ and $\Sigma_{1}^G=\Sigma_{2}^G=1.0428$ , while the SDGT2 calculation produced $\mT ^G=8.4\times10^{-22}$ with $\Sigma_{1}^G=\Sigma_{2}^G=0.99961$ . The CPU time required to transform this graph using our implementations of the SDGT1 and SDGT2 methods was and seconds, respectively.

To calculate the rates at temperatures in the interval reliably we used an implementation of the SDGT2 method compiled with quadruple precision (SDGT2Q) (note that the architecture is the same as in other benchmarks, i.e. with 32 bit wide registers). The points in the blue dataset in Figure 4.16 are the results from 4 SDGT2Q calculations at temperatures $T\in\{0.10,0.35,\dots,0.75\}$ .

**Figure:** Arrhenius plots for four digraphs of varying sizes (see Table 4.1) created from a sample of the PES for the LJ $_{55}$ cluster. is the rate constant corresponding to surface-to-centre migration of a tagged atom. Calculations were conducted at **$T\in \{0.3,0.35,\dots ,0.7\}$** using the SDGT method (green) and **$T\in \{0.1,0.15,\dots ,0.25\}$** using SDGT2Q (blue). For each of the digraphs the calculations yielded essentially identical results so data points for only one of them are shown.
$\begin{psfrags} \psfrag{0} [br][br]{$0$} \psfrag{2} [bc][bc]{$2$} \psfrag{... ...aphics[width=.5\textheight]{markov/arrhenius/LJ55/arrhenius.eps}} \end{psfrags}$

By using SDGT2Q we were also able to improve on the low-temperature results for LJ $_{38}$ presented in the previous section. The corresponding data is shown in blue in Figure 4.14.