Locally adaptive method to define coordination shell

The coordination shell of a particle represents its nearby particles. It is an important quantity to explain the structure of condensed-phase systems¹ and is necessary input for numerous order parameters.^2–6 The associated coordination number N_c is the number of particles in the shell. Despite their importance, there is no robust, general definition. The simplest approach is nearest-neighbor analysis whereby the coordination shell is assumed to comprise the nearest N_c particles. This has the obvious disadvantage, though, of requiring some existing knowledge of the system. Nonetheless, it proves useful and reasonably accurate in certain cases, such as approximating N_c by the value in an appropriate crystal structure.

The most commonly used method extracts N_c from the radial distribution function g(r), derived either from x-ray or neutron diffraction or computer simulation. This is done in a number of ways: placing a cut-off at the first minimum of g(r) or r²g(r)⁷ or symmetrizing the first peak of rg(r) or r²g(r).^7–11 Of these, the g(r) cut-off (GC) is the most popular, despite its often abrupt nature and possible lack of a minimum. Another drawback of any g(r)-based method, as for any parametrization, is that it is only applicable to the particles, conditions, and structures from which it is derived. In general, a different g(r) is required for every type of interaction, set of conditions, or local structures, a number which grows combinatorially with system complexity. Applying a given cut-off elsewhere is an unknown approximation. Even then, g(r) requires averaging over a sufficient number of local structures in order to produce a smooth curve with a well-defined minimum. Yet this very averaging means that any derived cut-off is mean-field and only approximately able to account for local variations in the structure typical of disordered systems.

The only viable way to avoid such problems is to eliminate mean-field parameters altogether and derive the structure in a locally adaptive way. The most longstanding and popular such method uses Voronoi polyhedra which are the regions of space closest to each particle. Particles are defined as neighbors if their polyhedra share a common surface.^1,12–15 However, at non-zero temperatures the resulting N_c typically averages ∼14,^14,15 as predicted by Bernal,¹³ regardless of whether the system is a crystal or liquid. This is clearly in excess of the expected N_c = 12 for the fcc crystal.¹⁵ A recent method with a system-independent parameter is the solid-angle nearest-neighbor (SANN) algorithm¹⁶ which places the cut-off distance at the point where the sum of the solid angles of nearest neighbors equals 4π. Equivalently, a particle lies beyond the cut-off if its distance exceeds the sum of distances to all closer neighbors divided by N_c − 2. SANN yields reasonable N_c for liquids but has a weak dependence on density. Other proposed algorithms are based on the occlusion of two particles by an intermediate particle. Two of these still require a pre-defined distance cut-off parameter^17,18 while a third defines two particles as being blocked if there is a triangle lying between the two particles and whose vertices are all shorter than the particles’ separating distance.¹⁹ 

Our interest in parameter-free methods to determine N_c arose in the case of hydrogen bonds in liquid water.²⁰ Taking a donor’s strongest acceptor except when there is a stronger, competing donor-acceptor interaction eliminates the cut-off parameters commonly used in hydrogen-bond definitions, except that a hydrogen’s N_c cannot exceed one. Even though criteria have been found to identify bifurcated hydrogens²¹ and N_c for aqueous cations,²² a more general approach is needed for any system. Our solution presented here is the relative angular distance (RAD) algorithm. The idea of RAD is that a particle lies within a coordination shell as long as no other particle blocks it and no closer particle is itself blocked. The ability of a particle to block another particle from lying in the central particle’s shell is assumed to depend on three quantities: the solid angles that the two particles subtend relative to the central particle and the cosine of the angle that they subtend at the central particle. RAD is applied to systems of Lennard-Jones argon atoms arranged in simple cubic (sc), body-centred cubic (bcc), face-centred cubic (fcc), hexagonal close-packed (hcp) crystals, the liquid over a range of temperatures and densities, and a representative gas. Comparisons with GC N_c values are used to assess RAD’s performance.

RAD determines N_c in the following way. With respect to particle i, particle j is considered unblocked if for all particles k

where θ_jik is the angle subtended at i by j and k as shown in Figure 1, Ω_ij is defined as the solid angle subtended by the cone whose apex is located at particle i and which inscribes a circle centered at particle j with diameter σ_j, and Ω_ik is defined similarly. Equation (1) may be simplified by deriving an expression from the equation for the solid angle

where α_ij is the angle subtended by particle j from the center of particle i. Making the small angle approximation sin(α_ij/2) ≈ (σ_j/2)/2r_ij, where σ_j is the diameter of particle j and r_ij is the distance between the centers of particles i and j, leads to the equation

Substituting this into Eq. (1) and cancelling the factors of 2 yields

which, for equal-sized particles, as in Figure 1, becomes

The RAD rule is that particle j is deemed to be in i’s shell if j and every particle closer to i are unblocked. RAD is efficiently and non-recursively accomplished by working outwards along the sorted list of particle i’s nearest neighbors j, comparing only closer particles as blockers k, and ending the search when a blocked particle is found.

Equation (5) alone without the requirement that closer particles be unblocked is termed RAD_open. To demonstrate the appropriateness of the squared power in three dimensions, the effect of varying r’s power is examined. Similar to the SANN algorithm,¹⁶ RAD is asymmetric in that two particles may disagree as to whether they are within each others’ shells. The effect of this is examined by two more methods: RAD_and, where a disagreement causes neither particle to be in the other’s shell, and RAD_or, where a disagreement causes both particles to be in each other’s shell. Concerning computational efficiency, because θ_jik depends on three particles, RAD would be O(N³) in computing performance. However, because only ∼N_c nearby particles need be considered for the coordination shell or as potential blockers, the usual nearest-neighbour O(NlogN) scaling is achievable.

All simulations were performed using the sander module of AMBER 14. Lennard-Jones parameters for argon²³ have σ = 3.404 Å and ϵ = 0.9960 kJ mol⁻¹. For the sc, bcc, fcc, and hcp crystals, a cubic lattice of 216 unit cells was created, and argon atoms were placed at suitable lattice points with the shortest spacing set to σ, giving 216, 432, 864, and 864 atoms for the sc, bcc, fcc, and hcp systems, respectively. Because Lennard-Jones particles do not form stable sc or bcc crystals, restraints of 2 kcal mol⁻¹ were applied in those cases. Systems were minimized with 500 steps of steepest descent, followed by equilibration for 40 ps using periodic boundary conditions under constant volume and constant temperature. Structures were then collected every 1 ps for 1 ns for analysis. All crystals were run at a temperature of 60 K, and liquids, melted from the fcc crystal, were run at 40 K intervals from 140 K to 300 K and at densities of 0.6 σ⁻³ and 0.8 σ⁻³. A gas-phase simulation was also run at a temperature of 300 K and density 0.012 σ⁻³. GC N_c values were calculated by integrating r²g(r) up to the cut-off at the first minimum of g(r). To further assess the transferability of g(r) cut-offs, the cut-offs derived for the liquid and fcc crystal were applied to all systems.

Table I compares simulation-derived N_c values for RAD and GC together with some method variants for the sc, bcc, and fcc crystals, liquid, and gas. The direct application of RAD to the sc, bcc, fcc, and hcp crystal structures gives N_c = 6, 14, 12, and 12, respectively. These are in line with expectations, although bcc has an alternative value of N_c = 8 if only the closest shell at distance σ is considered, which is marginally closer than the second shell at $( 2 / 3 ) σ$⁠. This distinction weakens upon the inclusion of thermal fluctuations and disorder. N_c values for the hcp crystal are indistinguishable from those for fcc. Figure 2 displays the N_c histograms of the GC and RAD methods, while Figure 3 illustrates each system’s g(r) and g(r) for only the RAD coordination shell. The respective g(r) cut-offs are 4.1, 3.7, 4.1, 4.1, 5.4, and 7.0 Å for the sc, bcc, fcc, hcp, liquid, and gas cases. Because there is no well-defined minimum for the gas, the cut-off is taken where g(r) reaches 1.

The RAD and GC N_c values for sc, fcc, and hcp agree closely, and their histograms have sharply defined peaks at N_c = 6 and 12. For bcc, the histograms are broader and there is a sizeable discrepancy between the methods, with N_c being 8.4 for GC and 13.4 for RAD. This occurs due to the two-shell nature of bcc. Taking the g(r) cut-off where g(r) drops to zero gives N_c = 14, in line with RAD. For the liquid, RAD gives N_c = 8.9, which is lower than the GC value of 12.6. The RAD g(r) maps well onto the first peak and smoothly tails off to zero at around the g(r) cut-off, explaining the difference in N_c between the two methods. N_c derived from the symmetrized r²g(r) first peak^7,9 is calculated to be 7.2, placing RAD’s N_c between those of the two g(r)-based methods. The RAD histogram rarely lies above 12 whereas the most probable GC value is 13. N_c using the 4.1 Å fcc cut-off (GC_fcc) and the 5.4 Å liquid cut-off (GC_liquid) starkly show the cut-offs’ lack of transferability, with the crystal’s N_c now 15.9–18.0 and the liquid’s only 5.3. The gas-phase RAD value of 2.8 is much larger than the GC value at 0.4. RAD variants with r⁻² replaced by r⁻¹, r⁻³, and r⁻⁶ demonstrate how blocking increases with power, that r⁻² is optimal with respect to the crystals, with r⁻³ being very similar. RAD_open is almost identical to RAD for the crystals, 0.6 higher for the liquid and 2.2 higher for the gas. Similarly, RAD_and and RAD_or match RAD for the crystals, ±0.4 for the liquids, and ±0.7 for the gas.

Figure 4 demonstrates the effect of temperature and density on N_c for the liquid. Both RAD and GC predict a sizeable reduction in N_c with density. RAD values steadily decrease with temperature while GC values are largely invariant, apart from fluctuations between two different cut-offs because of small variations in the position of the g(r)’s minimum.

RAD is a locally adaptive, predictive algorithm to classify structures because it can identify the coordinations of particles from a single configuration without any parameters or knowledge of the system. This contrasts with GC for which a suitable cut-off must first be derived for the conditions of interest from an ensemble of structures appropriate to the system being studied and sufficiently large to give a clear minimum at the required resolution. This cut-off is mean-field, not transferable, rather abrupt, and like any parameter has an associated error. RAD captures the variation of N_c in both density and temperature, and between the crystal and liquid, suggesting that it is a useful order parameter. In contrast, GC here only responds to density. The sensitivity of the cut-off at the minimum to temperature could be eliminated by further sampling but this only adds to the demands of GC. RAD predicts that the liquid has a lower N_c than the fcc crystal, argon’s stable crystalline phase, as opposed to GC which predicts a slightly higher value. A decrease is arguably more appropriate, given a liquid’s greater free volume and lower density. Furthermore, the RAD value lies between GC’s N_c and the lower N_c derived from symmetrizing the g(r)’s first peak.

The RAD inequality, Eq. (5), combines two intuitive and competing effects of proximity and co-linearity. The interpretation of the functions r⁻² and cosθ in inverse-square laws and the dot product, respectively, is long recognized. RAD may be readily generalized to particles with different sizes. A key feature of RAD compared to other methods is the exclusion of any particle if closer particles are blocked. This automatically identifies a local length scale for a coordination shell based on the shortest distance at which blocking occurs. Particles beyond this length scale cannot lie in the shell even if they are not blocked. Without this feature, as in RAD_open, the lowest possible value of N_c is 4, but with it N_c = 1. This difference is expected to be important in systems with widely varying density, such as the vapor-liquid interface, because distant particles are prevented from being in the shell of exposed particles at the surface. This is an advantage over SANN, which assumes a total solid angle of 4π for any particle, regardless of the local density, and thus overestimates N_c at the air-water interface.¹⁶ Table I makes clear that both r⁻¹ and RAD_open produce excessive N_c. RAD_and and RAD_or give very similar values to RAD and avoid the asymmetry whereby a particle lies in another’s shell but not the other way around. Asymmetry is expected to be more pronounced for particles of different sizes. The r⁻⁶ variant is based on the idea that blocking is proportional to the attractive pairwise Lennard-Jones force. Its low values of N_c, however, show that nearby particles have too much blocking power. The small-angle approximation in Eq. (3) performs well because the angle subtended by particles of equal size is only 30° at r = σ, the distance below which particles rarely approach. It may be better to remove this approximation in certain situations, such as for particles with very different sizes.

RAD is a general method to determine a particle’s coordination shell. Its application to collections of Lennard-Jones particles in crystalline, liquid, and gas phases at various temperatures and densities yields expected coordination numbers for the crystals and reasonable values for liquid and gas. Its variation with temperature and density makes it a suitable order parameter for tracking structural changes and phase transitions. While we have only tested the three main phases here, we anticipate that RAD should be general for all kinds of particles, conditions, structures, and dimensionalities.

This work was funded by BBSRC Grant No. BB/K001558/1.