An algorithm is presented to define a particle’s coordination shell for any collection of particles. It requires only the particles’ positions and no pre-existing knowledge or parameters beyond those already in the force field. A particle’s shell is taken to be all particles that are not blocked by any other particle and not further away than a blocked particle. Because blocking is based on two distances and an angle for triplets of particles, it is called the relative angular distance (RAD) algorithm. RAD is applied to Lennard-Jones particles in molecular dynamics simulations of crystalline, liquid, and gaseous phases at various temperatures and densities. RAD coordination shells agree well with those from a cut-off in the radial distribution function for the crystals and liquids and are slightly higher for the gas.
The coordination shell of a particle represents its nearby particles. It is an important quantity to explain the structure of condensed-phase systems1 and is necessary input for numerous order parameters.2–6 The associated coordination number Nc 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 Nc 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 Nc by the value in an appropriate crystal structure.
The most commonly used method extracts Nc 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 r2g(r)7 or symmetrizing the first peak of rg(r) or r2g(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 Nc typically averages ∼14,14,15 as predicted by Bernal,13 regardless of whether the system is a crystal or liquid. This is clearly in excess of the expected Nc = 12 for the fcc crystal.15 A recent method with a system-independent parameter is the solid-angle nearest-neighbor (SANN) algorithm16 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 Nc − 2. SANN yields reasonable Nc 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 parameter17,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.19
Our interest in parameter-free methods to determine Nc arose in the case of hydrogen bonds in liquid water.20 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 Nc cannot exceed one. Even though criteria have been found to identify bifurcated hydrogens21 and Nc for aqueous cations,22 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 Nc values are used to assess RAD’s performance.
A. RAD algorithm
RAD determines Nc 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)/2rij, where σj is the diameter of particle j and rij 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 RADopen. 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,16 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: RADand, where a disagreement causes neither particle to be in the other’s shell, and RADor, 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(N3) in computing performance. However, because only ∼Nc nearby particles need be considered for the coordination shell or as potential blockers, the usual nearest-neighbour O(NlogN) scaling is achievable.
B. Molecular dynamics simulations
All simulations were performed using the sander module of AMBER 14. Lennard-Jones parameters for argon23 have σ = 3.404 Å and ϵ = 0.9960 kJ mol−1. 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−1 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 σ−3 and 0.8 σ−3. A gas-phase simulation was also run at a temperature of 300 K and density 0.012 σ−3. GC Nc values were calculated by integrating r2g(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 Nc 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 Nc = 6, 14, 12, and 12, respectively. These are in line with expectations, although bcc has an alternative value of Nc = 8 if only the closest shell at distance σ is considered, which is marginally closer than the second shell at . This distinction weakens upon the inclusion of thermal fluctuations and disorder. Nc values for the hcp crystal are indistinguishable from those for fcc. Figure 2 displays the Nc 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.
|Method .||sc .||bcc .||fcc .||Liquid .||Gas .|
|Method .||sc .||bcc .||fcc .||Liquid .||Gas .|
The RAD and GC Nc values for sc, fcc, and hcp agree closely, and their histograms have sharply defined peaks at Nc = 6 and 12. For bcc, the histograms are broader and there is a sizeable discrepancy between the methods, with Nc 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 Nc = 14, in line with RAD. For the liquid, RAD gives Nc = 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 Nc between the two methods. Nc derived from the symmetrized r2g(r) first peak7,9 is calculated to be 7.2, placing RAD’s Nc between those of the two g(r)-based methods. The RAD histogram rarely lies above 12 whereas the most probable GC value is 13. Nc using the 4.1 Å fcc cut-off (GCfcc) and the 5.4 Å liquid cut-off (GCliquid) starkly show the cut-offs’ lack of transferability, with the crystal’s Nc 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−2 replaced by r−1, r−3, and r−6 demonstrate how blocking increases with power, that r−2 is optimal with respect to the crystals, with r−3 being very similar. RADopen is almost identical to RAD for the crystals, 0.6 higher for the liquid and 2.2 higher for the gas. Similarly, RADand and RADor 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 Nc for the liquid. Both RAD and GC predict a sizeable reduction in Nc 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 Nc 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 Nc 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 Nc and the lower Nc 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−2 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 RADopen, the lowest possible value of Nc is 4, but with it Nc = 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 Nc at the air-water interface.16 Table I makes clear that both r−1 and RADopen produce excessive Nc. RADand and RADor 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−6 variant is based on the idea that blocking is proportional to the attractive pairwise Lennard-Jones force. Its low values of Nc, 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.