Controlling disorder in two-dimensional networks

Two-dimensional networks are constructed by reference to a distribution of ring sizes and a parameter (α) which controls the preferred nearest-neighbour spatial correlations, and allows network topologies to be varied in a systematic manner. Our method efficiently utilizes the dual lattice and allows the range of physically-realisable configurations to be established and compared to networks observed for a wide range of real and model systems. Three different ring distributions are considered; a system containing five-, six- and seven-membered rings only (a proxy for amorphous graphene), the configuration proposed by Zachariasen in 1932, and a configuration observed experimentally for thin (near-2D) films of SiO2. The system energies are investigated as a function of the network topologies and the range of physically-realisable structures established and compared to known experimental results. The limits on the parameter α are discussed and compared to previous results. The evolution of the network structure as a function of topology is discussed in terms of the ring–ring pair distribution functions.

S Supplementary material for this article is available online (Some figures may appear in colour only in the online journal)

Letter
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. arrangements. In figure 1(a) similar ring sizes tightly cluster together, whereas in figure 1(c) they are dispersed throughout the structure. The rings comprising any given distribution may be arranged in a large number of configurations, only a subset of which may be physically-realisable, and so require an additional measure for characterisation. This is often provided by the Aboav-Weaire law, which may be written as where m n is the mean ring size about a central ring of size n.
The tendency of small rings (n < 6) to be adjacent to large rings (n > 6) is then captured in a single parameter, α (see, for example, the reviews of [17,21]). The relationship between α and the network topology can be seen again through figure 1, where increasing values of α accompany the increase in small-large ring pairings. There have been many experimental and computational studies which generate random networks and calculate the Aboav-Weaire parameter [6,[23][24][25]. These have found that the law holds well for a wide range of systems (with notable exceptions [26,27]). Despite large differences in the systems, many such networks display apparently similar properties with typical values of α in the range 0.16-0.30 [28]. Such studies are often confined to naturally occurring networks, and the exact physical meaning, and limits of, the Aboav-Weaire parameter remain unclear. The key to the present work is that methods are introduced in which both the ring distribution p n and Aboav-Weaire parameter are imposed on the network. The question being posed here is, therefore, given a pre-determined ring distribution (such as that obtained from experiment), how can these rings be arranged and what are the geometrical and physical limits? Three types of system are investigated. The first is a random network containing only five-, six-and seven-membered rings, a proxy for amorphous graphene (hereafter aG) and relatively well-defined. The second is based on the most famous example of a 2D network, Zachariasen's glass. Zachariasen was the first to represent systems such as SiO 2 as networks of percolating rings, and sketched a simple 2D configuration (reproduced in figure 1) purely in terms of the ring structure [22]. It is almost certain that Zachariasen paid little attention to p n and the detail of the nearest-neighbour ring distributions, his aim being to highlight how sketching bonds could define the network topology. Indeed, for his configuration p n is 'skewed' towards 5-rings and n = 5.90. The calculation of α is somewhat contrived for a small aperiodic sample such as his, but α = 0.31 if the six central rings are used and α = 0.34 if the edge rings are included, placing Zachariasen s glass at the far end of the experimentally-observed range of values. We therefore take a modified version of Zachariasen s distribution, where one 4-ring has been omitted from the original configuration to ensure n = 6, as an example of an 'extreme' distribution containing rings in the range n = 4-8. Finally we take an experimental distribution obtained from the study of silica bilayers which we anticipate as being more physically-realistic than than proposed by Zachariasen, with ring sizes in the range n = 4-10 [29].
We consider trivalent atomic networks, where each atom is bonded to exactly three others, as is the case for amorphous graphene. For SiO 2 the rings are constructed from Si-O-Si links, leading to wider observed ring distributions, but the oxygen linkages do not affect the ring topology. A network consists of N distinct ring sizes, which belong to the set R = {n 0 , n 1 , . . . , n N }, so that, for example, for aG (2) The atomic network does not, however, provide the most natural framework to consider the Aboav-Weaire law. This is because the atomic network is defined by the connectivities between neighbouring atoms, whereas the Aboav-Weaire law is concerned only with the connectivities between neighbouring rings. A more suitable approach is to consider the dual lattice, constructed by placing a node at the centre of each ring, and connecting the nodes of adjacent rings (highlighted for the hexagonal lattice in figure 2), such that the number of connections of the node is equal to the ring size. The dual lattice consists of tessellating triangles, and the atomic network can be reformed by placing an atom at the centre of each triangle.
To quantify the connectivity in the dual, we introduce the N × N matrix π, where each element, π nn , gives the proportion of rings of size n occupying adjacent positions to rings of size n. Each row of the matrix therefore corresponds to the individual ring distributions around a ring of size n, and is normalised so that n ∈R π nn = 1. The matrix π is not symmetric due to the differing number of connections in the rings n and n . They are however linked by the relationship: π nn = π n n n p n np n .
The mean ring size about each ring can be found using the summation which can be readily used in conjunction with the Aboav-Weaire law. Amorphous atomic networks are frequently generated using a 'bond-switching' algorithms in which Monte Carlo moves are accepted or rejected based on the geometry optimised energy of the system [30,31]. For this work, we have implemented a dual-space equivalent of this algorithm, where connections between nodes in the dual lattice are switched. Crucially however, moves are not accepted on the basis of the energy of the network, but rather its agreement with a target ring distribution and Aboav-Weaire parameter. Therefore, whereas bond-switching gives as its output the ring distribution and α parameter that correspond to physically reasonable structures, dual-switching takes these parameters as its input and evaluates the physical feasibility. This provides fine control over disorder in the network and allows for systematic investigation into the range of configurations available for a given system and their properties.
The dual-switching algorithm proceeds as follows. The starting structure is a periodic hexagonal lattice. For a proposed Monte Carlo move, a pair of edge sharing triangles are randomly selected, and the node connections transposed, incrementing the sizes of two rings and decrementing two others, as shown in figure 2. This move preserves both the total number of rings and n = 6. In two dimensions, this move is sufficient to access all possible network configurations [32].
An Aboav-Weaire fit is calculated from the dual to give values of α and μ. The network properties are compared to the target input values ( p t n , α t and µ t ), using the following cost function: where k is a scaling constant. In the cost function the relative difference must be used for the ring distribution, as the same accuracy is required for all p t n , which may differ by two orders of magnitude. This is not a concern for α t , which must also have the flexibility to take a zero value, and hence the absolute difference is used in the first term. The move is accepted with probability given by the Metropolis condition: where ∆χ is the difference in cost functions before and after the proposed move, and T is a 'temperature'. Moves are proposed until the network has converged to the target properties and the cost function is zero. The dual-switching algorithm generates a single network configuration with prescribed properties from a single calculation. It is in effect a minimisation routine, searching down the cost gradient to ultimately find network configurations which satisfy χ = 0. As is the case with local optimisation techniques, steps must be taken to avoid becoming trapped in local minima, resulting in non-convergence. This is achieved through selection of the parameters k and T. The parameter k changes the relative costs of satisfying the α t and p t n conditions, and must be chosen so that neither is overweighted (which would lead to deep local minima). The parameter T, an effective temperature, controls the proportion of 'uphill' moves which are accepted. Some temperature is required to overcome local minima, but if set too high the algorithm will no longer move consistently 'downhill' in cost and the search becomes effectively random-invariably leading to non-convergence. Values for k and T can be determined from a parameter search checking for convergence of target systems; but we found k = 10 and T ∼ 10 −4 to be appropriate for systems of the type and size described in this work.
Although there is no requirement for geometry optimisation within this algorithm, the energy of the system provides a measure of the physical feasibility of a given structure. Instead of continually generating and minimising the atomic network, which would be computationally expensive, we propose that a reasonable atomic configuration can be generated through optimisation of the dual lattice alone. Due to the triangular nature of the dual, a simple harmonic potential (as detailed in online supplementary section SI 3 (stacks.iop.org/ JPhysCM/30/50LT02/mmedia)) can be used for local optimisation after every accepted Monte Carlo move. After the Monte Carlo process has finished, the dual can be globally optimised before being converted to the atomic network and minimised with a more accurate atomic potential if necessary, which in this work was chosen as the Keating potential (see SI 3) [33]. All optimisations were performed at fixed density, although this could in principle be allowed to vary.
This method was applied first to aG networks. This system represents a useful framework for investigating the Aboav-Weaire law due to the presence of additional constraints which make them highly controllable. As a consequence of Euler's law the proportion of 5-and 7-rings must be equal, which gives a trivial relationship between the second moment and proportion of 6-rings, As stated above, the π matrix contains all the information needed to calculate α. For the aG system owing to the row normalisation, constraint (3) and by imposing the Aboav-Weaire law, it can be fully expressed in terms of only two elements, as shown in SI 1. These can be chosen to be π 55 and π 57 , which capture the tendency for the small rings to be adjacent to small and large rings respectively. It is also shown (SI 1) that, where ∆ 57 = π 57 − π 55 . Therefore, for the aG system, α is linearly related to the difference in 5-7 and 5-5 ring proportions, and hence can be described as a well-defined disordered system.
To find the range of accessible α values for the aG system, periodic networks containing 10 000 rings were generated with 0.1 p 6 0.9. The aim of these simulations was to try and probe the geometrical limits of α, and so a high number of Monte Carlo steps was used, 10 9 , without the need for geometry optimisation. Figure 3(a) shows the range of accessible α values as a function of p 6 . The upper limit, α max , appears a relatively weak function of p 6 whilst the lower limit, α min , shows a much stronger dependence. The range of accessible values, ∆α = α max − α min , broadly mirrors the system entropy (see SI 2), although there is deviation around p 6 = 1/3 (the presence of which is being investigated). Furthermore, α min are more negative than previous theoretical predictions which suggest limits of α min = −µ/36 [34] and α min = −µ/18 [35] respectively, the former being equivalent to the 'topological gas' [36] for which m n (equation (1)) is simply related to μ. Note that α does not necessarily tend to zero as p 6 −→ 1. For the case of a pair of 5-and 7-rings, α will tend to zero if the rings are 'isolated' (i.e. surrounded by 6-ring nearest-neighbours) but not if they are themselves nearest-neighbours. In some sense considering these problems in terms of a low concentration of defects is useful as the nearest-neighbour mean ring sizes can be expressed analytically, however the Aboav-Weaire law is only really appropriate when considering truly disordered networks.
To explore the structural properties of the aG networks at different values of { p 6 , α}, 100 periodic networks containing 10 000 rings, were constructed for p 6 = 0.2, 0.4, 0.6, 0.8. These simulations were performed with geometry optim isation and so also provide information on the physical limits on α. Snapshots of three example networks are given in SI 6. Figure 3(b) displays the mean and standard deviation of the total potential energy for each p 6 atomic network across a range of α. It can be seen that the energy minimum in each case is only weakly dependent on the value of p 6 , varying from The variation in the maximum and minimum values of α (denoted α max and α min respectively) for which networks can be generated, as a function of p 6 . The figure also shows the mixing entropy, ∆S/k B , as described in SI 2. (b) Energies for the aG system generated as a function of both α and p 6 (key as shown in the legend), the configurations derived from the original Zachariasen [22] and SiO 2 -bilayers [29]. α 0.23 at p 6 = 0.8 to α 0.27 at p 6 = 0.2, and close to the value of α seen across many natural systems. Whilst there is little cost for small deviations from the minimum, decreasing α rapidly incurs a relatively large energetic penalty. Figure S1 shows the analogous energies when minimising through the dual lattice alone. The curves have a very similar form with the minima aligned, suggesting that working in dual-space can be sufficient to capture all system properties, with a much lower computational overhead.
To further quantify any ordering imposed on the generated configurations radial distribution functions, g nn (r), are constructed which are coloured by ring size, {n, n } and provide information regarding the mean average distances between the centres of rings of different sizes. Figures 4(a) and (b) show the rdfs for the 5-5 and 5-7 ring pairings, g 55 (r) and g 57 (r) respectively (the remaining functions are shown in the supplementary information, figure S2). As is consistent with the intuitive meaning of α, increasing α causes a reduction in intensity in the first peak of g 55 (r) and a concomitant increase in intensity in the first peak of g 57 (r), as 5-5 adjacencies are replaced with 5-7. In addition, the position of the first peak shifts to smaller r as α is reduced, reflecting both the increased distortion in the rings and the deviation from the ideal 120° bond angle, which translate to the higher observed potential energy. Figures 4(a) and (b) also show significant structural evolution beyond the nearest-neighbour length-scale. As α becomes more positive peaks emerge in g 55 (r) at r/r • 55 ∼ 1.8 and ∼2.3. An increase in α corresponds to a greater tendency for 7-rings to be nearneighbours to 5-rings and, in turn, increases the probability of the same 7-ring having a second 5-ring near-neighbour. In simple geometric terms the second 5-ring can occupy three possible sites around the 7-ring (see figure S3), the non-adjacent positions corresponding to the developing peaks. Note that one might naively assume that driving α to more positive values would tend to eliminate the nearest-neighbour 5-5 spatial correlations. However, figure 4(a) indicates this not to be the case, reflecting the balance between retaining these units and facilitating nearest-neighbour 5-7 ring interactions via the formation of 5-7-5 triplets.
Similar analysis was performed on 100 generated Zachariasen and SiO 2 networks, examples of which can be found in SI 7 and SI 8 respectively. We note here that although our algorithm requires the fit to equation (1) to be exactly linear for the aG system, for broader ring distributions this is no longer the case. However, for the Zachariasen configuration the linear regression (R 2 ) coefficient was always in excess of 0.995, and for the silica the average R 2 was 0.979, representing a very good fit. Figure 3 shows the energies of both the Zachariasen and SiO 2 systems as a function of α. Both cases resemble those for the aG with energy minima at α ∼ 0.25. The silica curve shows smaller curvature reflecting the broader distribution of ring sizes whilst the Zachariasen curve shows a greater curvature reflecting the 'extreme' (i.e. physically unrealistic) nature of the distribution. In addition it proved difficult to generate low α configurations (α < −0.1) for the Zachariasen network. Figures 4(c) and (d) show two key rdfs for the Zachariasen configuration, g 44 (r) and g 88 (r), highlighting the spatial correlations between the smallest and largest rings in the system. The effects of changing α on g 44 (r) are dramatic with strong nearest-neighbour clustering at negative values. In this case, however, the nearest-neighbour 4-4 correlations do vanish at high α as 4-8 nearest-neighbour correlations dominate but the 8-ring is large enough to accommodate up to four 4-ring nearest-neighbours without any 4-4 neighbouring pairs. Again this is demonstrated through the next nearest neighbours by the 8-4-8 peak developing at ∼1.4.
The work discussed here has focused exclusively on 2D networks with one nearest-neighbour length-scale, corresponding to a monatomic system. Future possible developments vary in complexity. For example, introducing a variety of atom identities (and hence nearest-neighbour lengths-scales) is effectively trivial if trivalency is maintained. Introducing a dispersion in valency is more involved, as it removes the simple relationship between the atomic network and dual lattice. However, less computationally-efficient methods which do not require the dual lattice could be employed, or the atomicdual relationship further investigated. Similar comments apply to extending the method to three-dimensions, with the addition that work would be required in effectively developing a three-dimensional version of the Aboav-Weaire law (see, for example, [37]) Our key objective was to highlight the development of methods to generate configurations in which both the distribution of ring sizes and their topological arrangements can be specified. We have presented a method to generate amorphous structures in which the ring distribution and nearestneighbour correlations are finely controlled. The method is robust, as demonstrated by the construction a wide range of configurations for three very different example ring systems, and is consistent with experiment in the determination of the Aboav-Weaire α parameter for optimal network stability. Furthermore, the meaning of the α parameter has been quantified in terms of elements in a matrix, π, and its effect on the short and medium range order of networks illustrated through the calculation and analysis of partial radial distribution functions in the dual lattice.