A generalized scheme for characterizing orientational correlations in condensed phases of high symmetry molecules: SF$_6$ and C$_{60}$

The orientational correlation scheme introduced earlier for tetrahedral molecules is extended for being able to classify orientational correlations between pairs of high symmetry molecules. While in the original algorithm a given orientation of a pair of tetrahedral molecules is characterized unambiguously by the number of ligand atoms that can be found between two planes that contain each centre and perpendicular to the centre-centre connecting line, in the generalized algorithm, the planes are replaced by cones, whose apex angles are set according to the symmetry of each molecule. To demonstrate the applicability of the method, the octahedral-shaped SF$_6$ molecule is studied in a wide range of phases (gaseous, supercritical fluid, liquid and plastic crystalline) using classical molecular dynamics. By analyzing the orientational correlations, a close-contact region in the first coordination shell and a medium-range order behaviour are identified in the non-crystalline phases. While the former is invariant to changes of the density, the latter showed longer-ranged correlations as density is raised. In the plastic crystalline state, fluorine atoms are oriented along the lattice directions with higher probability. To test the method for icosahedral symmetries, the crystalline structures of room temperature C$_{60}$ is generated by three sets of potentials that produce different local arrangements. The novel analysis provided quantitative result on preferred arrangements.


Introduction
A most important issue in the structure analysis of disordered materials is to find suitable and simple representation(s) of the studied material that can emphasize its characteristic property. Considering orientational correlations in condensed phases of high-symmetry molecules (including the ones whose shape belong to the most symmetrical ones called 'Platonic solids': tetrahedron, octahedron, cube, icosahedron and dodecahedron), there is a lack of suitable descriptions: atomic radial distribution function (RDF, g(r)) analysis is a frequently applied method, but it does not take into account the symmetry of molecules in a direct way. On the other hand, the application of symmetry adapted functions is restricted mainly to crystalline phases. In 2007 Rey [1] introduced a scheme that classifies any mutual arrangement between pairs of tetrahedral molecules to an ideal arrangement. In order to characterize a given arrangement, two parallel planes are taken, each contains one of the molecular centres and perpendicular to the line connecting them. Then, the classification is performed by the number of corner (or ligand) atoms of each molecule are between the two planes. This way, six different kinds of ideal arrangements have been defined, covering all possible orientations, as corner to corner(1:1), corner to edge(1:2), corner to face(1:3), edge to edge(2:2), edge to face(2:3) and face to face (3:3). Thus, the probability of each arrangement can be determined for any centrecentre distances. The advantage of this method that it provides quantitative results and represents correlations continuously in distance rather than providing an average over a coordination shell. The method has contributed to the elucidation of the structure of liquids consisting of regular [2][3][4][5][6] and distorted tetrahedral shapes of molecules [7,8] and recently, of ionic conductors [9].
Until now, this method has not been used for pairs of molecules whose shapes are non-tetrahedral Platonic. Note using the parallel planes as characterizing objects would provide useless results: only one class is formed as the non-tetrahedral Platonic solids possess inversion symmetry, which assigns exactly half of the corner atoms oriented to the other molecule in almost any case. Thus, the method should be generalized in another way: by keeping each class (corner, edge and face) and introducing orientation cones instead of planes as discussed in Section 2. To demonstrate the applicability of this approach, orientational correlation analysis in condensed phases of the octahedral-shaped sulfur-hexafluoride (SF 6 ) and room temperature crystalline phase of the truncated icosahedral-shaped C 60 (based on various potential models) are selected.
Sulfur hexafluoride (SF 6 ) is regularly used as an insulator in high-voltage circuit breakers [10]. Apart from the technical importance, the high symmetry of the molecules and a peculiar body-centered plastic crystalline phase attracted much attention from a scientific point of view [11][12][13][14][15][16]. This phase is stable between the point of sublimation at 223K and phase transition to the lowtemperature monoclinic phase [17] occurs at 96K under ambient pressure [12]. While the centre of mass of the molecule possesses translational symmetry, apparent molecular rotations are only hindered by neighbouring molecules. Neutron diffraction [11,12,17] and molecular dynamics (MD) modelling [14,15] stud-ies usually discussed correlations relative to the ordered crystalline lattice and concluded that fluorine atoms oriented close to the cubic cell axis. While the orientation of the molecule relative to the lattice is frequently taken into account in MD studies addressing the low-temperature phase [18] and the ordered-plastic phase transitions [19,20], there are only a few MD [14,15], and one [16] Reverse Monte Carlo (RMC) studies that consider local atomic arrangements in terms of atomic RDF-s. RMC was developed originally for liquids and amorphous materials [21] and provide sets of atomic configurations consistent with measured diffraction data without the need to use any potential.
The gaseous, liquid and supercritical fluid states have triggered only a few diffraction studies [22,23]. MD simulations are compared to experimental properties of these phases by Dellis and Samios [24], using potential sets optimized for describing supercritical fluid [23] and crystalline phases [13,25]. Strauss et al. [23] used RMC modelling, and their description was focused on atomic RDFs. According to my best knowledge, there has not been any investigation performed yet that would take into account the symmetry of the SF 6 molecule and aim at mutual orientations both in liquid and crystalline phases.
C 60 (fullerene, buckyball) has numerous phases depending on temperature and pressure [26,27]. At low-temperature and ambient pressure, it forms a simple cubic molecular crystal with some disorder. At 260K molecules start to rotate, thus inducing a phase transition [28] to the orientationally disordered face-centered cubic structure, stable at and above room temperature. Early diffraction studies [29][30][31] clarified the structure of the molecule: they observed a difference between single and double bonds ('D') separating pentagon ('P')-hexagon ('H') and hexagon-hexagon faces, respectively. Concerning the room temperature phase, powder diffraction and RDF (as derived from powder diffraction) studies suggested free and independent rotation [28] and a lack of intermolecular ordering [30]. On the other hand, single-crystal studies provided evidence of a slightly (within 20% intensity modulation) preferred orientation of molecules [32] with more pentagonal than hexagonal faces turning towards nearest neighbour directions [33], and of local ordering of molecules by diffuse scattering [34][35][36][37].
To describe the structure and take into account orientational ordering, various potential sets have been developed, as reviewed by Launois et. al [38,39] and Chaplot et. al [40]. These potential sets spread from simple van der Waals interactions on carbons [41,42], additional van der Waals centres on doublebonds with charges [43], additional charges on single and double-bonds [44], split charge on double-bonds [36,45], anisotropic van der Waals [46], mean-field [47] to ab-initio [48]. The main reason to overcome the simple van der Waals potential is that it does not produce the experimentally observed low-temperature phase [49]. Among the reviewed forcefields, none of them is able to describe all of the observed properties [39]. Orientational correlations, for instance, depend significantly on the applied forcefield [40] Details of the MD simulations and calculation of the orientational correlations are presented in Section 3. Section 4 contains the multiphase (gaseous, liquid, fluid, plastic crystalline) investigation of orientational correlations of oc-tahedral SF 6 using the extended method. It is followed by the quantitative analysis of short-range order by different forcefields of room temperature C 60 , consisting of truncated icosahedral molecules. Finally Section 5 summarizes the results.

Generalization of the classification scheme
To extend the classification scheme introduced by Rey[1] originally for tetrahedral molecules, the question of how many corner atoms of each molecule is between two planes perpendicular to the connecting line of the molecular centres should be replaced. Instead of perpendicular planes, we should consider orientational cones with proper angle constructed by supposing: 1. The directions of corner atoms from the centre are distributed evenly on the surface of a sphere. 2. The average number of corner atoms from one molecule inside the cone should be 2 when the molecule is oriented randomly. 3. Only 1, 2, or 3 corner atoms are allowed to be inside the cone.
The first two conditions determine the half apex angle (γ) of the cone that belong to the molecule possessing N corner corner atoms. According to them, the ratio of the solid angle of the orientation cone and the total solid angle should be equal to the ratio of the average number of corner atoms inside the cone at random orinetation of the molecule and N corner : after some re-arrangement: is obtained. This corresponds to the definition of a plane in the case of regular or slightly distorted tetrahedral molecules. The next task is the determination of the random orientation probabilities for each situation: the third condition limits them to three cases only. However, the second condition fixes the average to be 2, so the P i probabilities of having i number of corner atoms inside a randomly directed cone for one and three atoms should be equal. Thus, only a single P i probability parameter exactly determines the other two: The derivation and calculation method of P i with the help of corner-centrecorner bond angles is presented in SI 1.
Combining these probabilities for pairs of molecules according to the appendix of [1], the asymptotic values of each of the 6 groups can be obtained for molecules having tetrahedral, octahedral and icosahedral symmetries (see TABLE 1). Table 1: Asymptotic probabilities of each groups for different symmetry of molecules obtained by taking into account the p 1:1 = p 3:3 = P 2 1 , p 1:2 = p 2:3 = 2P 1 P 2 , p 2:2 = P 2 2 and p 1:3 = 2P 2 1 relations for single component molecular pairs (according to [1]). P 1 is the probability to find one corner atom inside the cone, which is determined on the basis of angular separation of the number of corner atoms (Ncorner) and half apex angle of the cone (γ).
molecule N corner cos γ γ P 1 = P 3 P 2 = 1 − 2P 1 The discussion of the studied system using this formalism is facilitated by the fact that each of the 6 ideal arrangement is broadly used (for instance octahedral molecules the figure 1 of REF [14] contains almost all ideal arrangement except the edge to corner) and a given arrangement of a molecular pair is categorized unambiguously into one of the 6 idealized ones.

Simulations
Classical MD simulations with flexible molecules have been performed by the GROMACS [50] package using version 2016.3. All simulations used Lennard-Jones interaction parameters and partial charges to model intermolecular and distant intramolecular dispersive and Coulombic interactions. Except for the case of one forcefield where the molecules are kept rigid, intramolecular interactions between adjacent atoms are modelled with harmonic bond stretching, angle bending and improper dihedral potentials.
Concerning the case of room temperature phase of fullerene, three potential sets have been used to model intermolecular interactions: the 'Girifalco' [41] forcefield with Lennard-Jones potential on each carbon atom; the 'Lu' [44] forcefield with modified Lennard-Jones parameters on the atoms and additional charges on single and double-bonds; and the 'Sprik' [43] forcefield with charges and Lennard-Jones potentials on each carbon atom and double-bonds, with Lorentz-Berthelot combination rules. According to 'Sprik', molecules were kept rigid, while in the former two forcefields the flexible potential set of Monticelli et al. [42] has been used.
Further details are presented in SI 2.

Calculating orientational correlations
In the saved trajectories of SF 6 , positions of sulfur and fluorine served as the central and corner atoms, respectively. For the C 60 particle configurations the icosahedral symmetry was reconstructed for each molecule. One of the possible generation of the C 60 molecule from a perfect icosahedron is to truncate each edge at a given ratio around each corner. As in a corner of the icosahedron 5 edges join, a pentagonal face is formed after the truncation. On the other hand, each triangular face of the icosahedron is truncated and a hexagonal face is created. Then, the non-truncated edges of the icosahedron are the double bond s of C 60 where two hexagonal faces join. As in the following, non-endohedral fullerenes are taken, the centres of them were calculated as the mean of atomic coordinates and the centre of each pentagon served as the virtual corner atom of the icosahedron at calculations of centre-centre RDF and orientations.
Following the symmetry reconstruction, the categorization of the orientations of each pair of molecules into the 6 categories as a function of centre-centre distances is performed using custom-written software. During the examination of the pairs, the numbers of corner atoms inside the cones centred on the centrecentre connecting line and having half apex angles of 70.53 o for SF 6 and 48.19 o for icosahedral representant of C 60 have been taken into account, according to the requirement of section 2. Next, each contribution in a given distance bin is normalized by the sum of them. Finally, the probabilities calculated from each of the 101 configurations are averaged.

Orientational correlations in condensed phases of sulfur hexafluoride
Neutron diffraction measurements at various thermodynamical conditions (published for gaseous [22]: 293.15 K and 10.2 bar, ρ=0.0021 Å −3 ; supercritical fluid [23]: 398 K; at 128 bar, ρ=0.0245 Å −3 ; at 155 bar, ρ=0.0289 Å −3 ; at 394 bar, ρ=0.0404 Å −3 ; at 1827 bar, ρ=0.0534 Å −3 and plastic crystalline [16,51] states: 190 K and 1 bar, ρ=0.0685 Å −3 ) allow us to observe the changes of orientational correlations in the framework of the modified algorithm. Similarly, a large variety of forcefields are available for this material. In SI. 4., these forcefields are tested if they are able to reproduce the available total scattering neutron diffraction data (calculation method of them are presented in SI. 3.) and density. Concerning these quantities, the 'Dellis' ('7sites' in [24]) forcefield of Dellis and Samios [24] was proven to be the best and provided a good agreement with experimental results. Note that this forcefield had been optimized originally against densities in the liquid and supercritical fluid regions, it also performed well close to the triple point (225 K and 10.1 bar, ρ=0.0532 Å −3 ) and in the plastic crystal phase that had not been parts of the original range of validation. Thus, in the following, results of MD simulations from this potential set are used.
FIGURES 1 and 2 contain the probabilities of each group with the molecular  centre-centre radial distribution function in 6 non-crystalline states: one gaseous (FIGURE 1a), 4 supercritical fluid (FIGURE 1a-e) and a close-to-triple point liquid (FIGURE 1f).
In order to obtain a general picture, we start with the examination of the gaseous state (FIGURE 1a). In this phase the colliding molecules determine the behaviour of the system, so we expect that the direct contact of atoms characterizes the short-lived orientational correlations. The molecular centre-centre RDF reflects this feature as only a broad first maximum is present with a maximum at 5.3 Å, which slowly decays at higher distances. Focusing on the orientational correlations, they definitely depend on the molecular centre-centre distance: the face-face (3:3) arrangement is frequent at the shortest distances, while farther away, the probability drops until it reaches its equilibrium value at 6 Å, slightly over the first maximum of the centre-centre RDF. Next, the edge-face (2:3) and further, the edge-edge (2:2) arrangement probabilities rise and have definite maxima at distances shorter than the maximum of the centre-centre RDF. The remaining corner-face (1:3), corner-edge (1:2) and corner-corner (1:1) arrangements have only the rising part in probability without a characteristic maximum: they reach their asymptotic values roughly at 6Å.
The observed behaviour has a simple geometric reason: the centre to a centre of a face, to mid of edge and to the corner atom distances are increasing, respectively. If we wish to pack atoms of molecular pairs as close as possible, the order of appearance of the arrangements should be 3:3, 3:2, 2:2 and 1:3, 1:2 and 1:1 as the centre-centre distance increases. A similar tendency is valid in most of the tetrahedral molecules [2] when this method is applied. However, in the present case, there are apparent differences from it, namely the intensity of 3:3 is not 1 at the shortest distances, but it is somewhat comparable to the probability of 2:3 arrangements. The explaination may be that there are only a few pairs in this region.
Extending the discussion to the supercritical fluid and liquid states (FIGURE 1a-f), we can conclude that up to about 5.5 Å orientational correlations do not alter in comparison with the gaseous state. This common behaviour is the direct consequence of the close contact of corner atoms as discussed above.
While the gaseous state does not show any significant correlations, neither in the centre-centre, nor orientations beyond the first coordination shell, the centre-centre RDF becomes more structured as the density of the system increases. Oscillations also appear (FIGURES 1 and 2) in the orientational correlations above 5.5..6 Å, which is slightly less than the upper limit of the first coordination shell. These oscillations are observed up to the upper limit of the second coordination shell in the two lowest density fluid states, up to the third maximum of centre-centre RDF for 394 bar and up to the upper limit of the fourth shell for the highest density states. It is important to note here that the simulated molecules were flexible and atomic coordinates are used without any further alignment during the calculation of orientational correlations, which slightly biased the asymptotic probability value of some groups (the largest deviation is about 0.002 for the 2:2 pairs).
Concerning the different kind of arrangement at medium-range, they do not show any preference for a specific correlation, but seemingly the probabilities of the 3:3 and the 3:2 arrangements are in opposite phase with probabilities of the 1:1 and the 1:2. Importantly, except the first common node point of probabilities of the 3:3 with the 1:1, of the 2:3 with the 1:2 and where the centre-centre RDF passing through unity, they coincide with each other. Moreover, the probability of arrangements having more corner atoms inside the orientational cones than the average 4 (the 3:3 and the 3:2) are in phase with, while arrangements indicate fewer corner atoms are in opposite phase with the centre-centre RDF. This observed medium-range order suggests a relationship with the tetrahedralshaped CCl 4 liquid[1], where similar behaviour was found. However, due to the remarkable short-range order correlations, the medium-range order starts there about from the third coordination shell and correlations are slightly shifted in phase with the centre-centre RDF (see FIGURES 4 and 5 of [1]). Thus, the two systems behave differently.
One of the advantages of the formalism just introduced is that it can be applied to not just disordered, but crystalline phases, as well. In the bodycentered cubic plastic crystalline state at 190 K, each sulfur atom has a definite equilibrium position, while fluorine atomic positions are disordered. This phenomenon implies that although equilibrium positions of sulfur atoms keep the translational symmetry, the molecular centre-centre RDF (FIGURE 3) expresses noticeable thermal displacements. Thus, the nearest neighbours in the <111> (or half of cubic diagonal), and second neighbours in the <100> (or side of the cube) directions of the crystal are overlapping in the radial representation, as their ideal values are close: 5.10 Å and 5.89 Å, respectively. Although the present model is probably more disordered than in reality (as discussed in SI. 4.), the overlapping behaviour was also observed earlier by Dove and Pawley [14] using another forcefield. Still, orientational correlations reveal some further features, allowing to separate the contributions of overlapping neighbours: in the <111> direction the 3:3 and 2:3 orientations dominate, meanwhile the 1:3 behaves like in the liquid states and correlations of 2:2 are less abundant. It is interesting to note that the probability of the 3:3 arrangements for close contact is high, while the 2:3 probability is more competitive with 3:3 for the studied liquids and gaseous states. At larger distances, the probabilities of 3:3 and 2:3 drop to the level of random orientations or even less, which shows the approximate limit between the first and second neighbours around 5.5 Å. This distance was the rough upper limit of the effect of close contact in the disordered phases. Nevertheless, the most significant characteristic of neighbours in the <100> directions is the broad maxima of 1:1 and 1:2 arrangements. In the <110> (or diagonal of a face of the cube) direction (between 7 and 9 Å) the 2:2 and 2:3 arrangements have maxima. Combining the arrangements found like the maxima of 1:1 to <100> and 3:3 or 2:3 to the <111> directions, the preferred position of fluorine atoms with respect to the lattice can be derived. For the first neighbour, 3 fluorines of each molecule should be placed inside the orientational cones around the cubic body diagonals, while at the same time only one fluorine is allowed to be inside the cone around the side of the cubic cell. the side of the cubes in the ideal case) cover directions of the three cubic sides, the fluorine atoms should be placed approximately toward the <100> direction from the sulfur positions with higher probability, which agrees well with earlier suggestions [12,15,16].
4.2. Short-range order orientation correlations in models of room temperature C 60 As a next step, the orientational correlation scheme is applied to quantitatively describe orientational correlations of neighbouring molecules in C 60 , as represented by three forcefields ('Girifalco' [41], 'Lu' [44] and 'Sprik' [43]) that had been proven to produce different kinds of short-range order.
First, we start the examination of the molecular centre-centre RDF (FIG-URE 4). Although some models show a larger centre of mass motion, the  most important feature of these RDFs, that the first and second neighbours are distinct. This allows separating orientational correlations unambiguously belonging to the different spheres, in contrast to the plastic crystal phase of SF 6 , where first and second shells overlap. Touching the broadening issue, it is the result of a different kind of modelling of the intramolecular structure. It was flexible for the 'Girifalco' and 'Lu' models, so the motion of each individual carbon atom is successfully equilibrated by the internal freedom, making the centre less intact from the changes. On the other hand, the 'Sprik' model was rigid, which resulted in a larger displacement of the centres from the equilibrium positions.
Focusing on the orientational correlations, the number of virtual corner atoms in the coordination cone should be assigned to the orientation of C 60 molecule. One corner atom (the "corner") in the coordination cone means the molecule is oriented by a pentagonal face ('P'), while two atoms ("edge") means between the centres of pentagons: the 'double-bond' ('D'). Finally, the "face" coordination realizes an orientation through the hexagonal faces ('H') of the molecule.
As earlier investigations focused on the short-range order behaviour, here also the first and second neighbours are examined. In order to quantify orientational correlations within a specific coordination shell, the probabilities of each group have been averaged over the given shell. The averaging has been made between the two limiting distances where g cc (r) intensity is 1/10 of the maximum intensity within the corresponding shell. Then the shell-averaged probabilities are divided by the corresponding random orientation probability (TABLE 1) for icosahedral molecules: Thus, its value over 100% indicates a preferred, while lower refers to a disliked orientation in TABLE 2. However, the probability of some arrangements  TABLE 1 are discussed together with tendencies on FIGURE 4. In the first coordination shell, which corresponds to the closest atoms to <110> directions in the crystalline system, the 'Girifalco' model shows a strong preference of 'HH', 'HD' and a slight preference of 'DD' orientational groups. This suggests that local arrangements try to avoid turning into any of the pentagonal faces of them and prefers hexagonal face -hexagonal face or doublebonds at closer and double-bond-double-bond orientation at distant centrecentre distances, which gives a hint about the rotation mechanism. Looking at the second coordination shell (<100> direction), the opposite 'PP' and 'PD' are favoured strongly by increasing centre-centre distances and the 'DD' slightly, but mainly at closer distances. Even though the preference of 'HD' and unlike of 'PP' and 'PD' found in accordance with earlier study [40], it contradicts to singlecrystal investigation [33] concerning the preference of pentagonal over hexagonal faces for nearest neighbours.
While the pentagonal coordination was disfavored in the first coordination shell by the 'Girifalco' model, the result by the 'Lu' forcefield shows a clear preference of it, agreeing with [33]: both 'PH' (mainly closer distances) and 'PD' (mainly higher distances) kinds of correlations are favoured and the 'PH' arrangement at short distance is more probable. The second shell shows less varied arrangements: only the 'DD' arrangement is preferred. This result is in accord with earlier investigation [40] about the first shell of 'PD' liked and 'PP' disliked.
Interestingly, the 'Sprik' forcefield shows a somewhat opposite kind of preference in comparison with the 'Girifalco', but the deviations from random values are moderate. The pentagonal face related correlations are preferred in the first shell in order of decreasing probability as 'PP', 'PD', 'PH', while in the second shell the hexagonal face related to 'HH', 'HD' and 'DD' correlations are the preferred ones. The 'PP' and 'PD' preference of this model is in agreement with an ab-initio study [48].

Conclusions
Based on the results the following statements can be made: (i) The orientational correlation scheme introduced originally by Rey [1] is extended to be capable of describing correlations between high symmetry molecules. The orientations of pairs of molecules are classified by the number of corner (aka ligand) atoms for each part found within orientational cones. The axes of the orientation cones correspond to the centre-centre connecting line with the apex at the centres. The extension is made by choosing its apex angle so that the number of corner atoms spans only from 1 to 3 and the average number over random orientation should be 2. Supposing even distribution of atoms, calculation methods for the half apex angle of the cones and asymptotic probabilities of each class are provided. (ii) In order to demonstrate the applicability of this method in various condensed phases, the octahedral-shaped SF 6 molecule is selected in its gaseous, supercritical fluid, triple point liquid and plastic crystalline phase. The modelling performed by the '7Sites' forcefield [24] provided the best agreement with density and total scattering diffraction data.
(iii) The detailed analysis of different orientational groups for SF 6 revealed a close contact region in the first coordination shell. It shows a specific orientation pattern of orientational arrangement by distance, similarly to that found for tetrahedral molecules. For non-crystalline states, a mediumrange ordered region is found that gradually extends with increasing density, even up to the 4th shell, as observed in high-density states. Its main characteristic feature is that the molecular centre-centre RDF is in phase with the probabilities of arrangements having more corner atoms within the orientational cones than the average 4. Although medium-range order appears for the tetrahedral CCl 4 too, it is different from the one observed here. In the plastic crystalline state, close contacts of 3:3, 2:3 correlations are found in the <111> directions, while the second neighbour in <100> showed the preference of 1:1 and 1:2 arrangements. This agrees with earlier studies, that fluorine atoms preferentially oriented in the <100> direction. (iv) For the room temperature crystalline C 60 , three models ('Girifalco' [41], 'Lu'[44] and 'Sprik' [43]) were tested in the framework of icosahedral symmetry, providing different short-range order. The probability of each arrangement is averaged over the first and the second coordination shells and divided by the corresponding random probability. Thus a quantitative description of orientational correlations is provided for each model, where agreement found in each case with earlier studies.
It is important to note that this method might be used not only to describe mutual orientational correlations between molecules belonging to the 'Platonic' solids (the un-mentioned cube and dodecahedron are the corner-face duals of octahedron and icosahedron, respectively) but also between their distorted forms or any object, as long as the three conditions in section 2 are fulfilled (like for trigonal planar or bipyramidal molecules, where N corner is 3 and 5, respectively). Also, it is possible to use the new algorithm for their mixtures (like tetrahedral and octahedral).

Acknowledgement
The author would like to thank the help of László Pusztai for discussions. The author is grateful that this project was supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences and the National Research, Development and Innovation Office(NKFIH), under grant no. KH130425.

Calculating the random P i probabilities
In order to determine the P i values for a given molecule, we should consider that each corner atoms are spaced more or less regularly. Thus each P i probability corresponds to the surface ratio on a unit sphere, where the axis might pass through and there are i number of atoms inside the cone (according to EQ. A2 of REF [1]). Now, let's fix a molecule and create a cone with γ half apex angle around each centre to corner atom axis. By doing this, we get a picture like in fig S1 around three corner atoms. If a molecular centre-centre connecting line passes Figure S1: Surface mapping on the unit sphere. The blue coloured arrows point toward ligand (corner) atoms from the centre of the molecule serving as axes of cones, which intersects the surface at lines marked by red circles. The light, middle and dark purple regions are surface areas containing 1, 2 and 3 ligands if the molecular centre-centre connecting line pass through them. To calculate surface ratios easy, individual surfaces of corner (p cI ; inside red circle), of edge (p eI ; surface limited by green lines) and of face (p f I ; surface limited by black lines) are introduced.
through one of these regions on a unit spheres, which might belong to 3, 2 or only 1 such cone(s), it will find the corresponding number of corner atoms inside the orientation cone. Thus, summarizing all of those individual surface elements, which belong to i number of cones, we get the corresponding probability.
For practical reasons, we introduce the corner (p c : surface ratio of a cone), the edge (p e : surface ratio of 2 overlapping cones) and the face (p f : surface ratio of three overlapping cones) surface ratios composed by the sum of each individual ratios (p x = p xI ) on the sphere (see fig S1). It is important to note, that p c contains areas, which belong to not just one cone and similarly p e contain areas, which belong to 3 cones. As a result, p c and p e are more than one.
Connecting these ratios to the P i probabilities, the P 3 is exactly equal to p f . On the other hand, the surface of P 3 is considered 3 times in the calculations of each p c and p e . Similarly, the surface of P 2 is also considered twice in p c calculation: Using EQ. 3 we get p c = 2 (S4) EQ. S4 is just the reformulation of EQ. 1, that the orientation cone should contain 2 corner atoms. The least 2 equation states, that either p f or p e should be calculated to determine all P i probabilities. Although, appendix of [1] calculated the p f I for tetrahedral molecules, in the following the p eI is calculated. It takes into account only pair of corner atoms, whose associated parameter (the cornercentre-corner bond angle) is more convenient to use in the case of non-regular shaped objects. In order to calculate the surface ratio, we should consider two cones with a common apex, both having the same half apex angle (γ) and their axes are separated by α I degrees from each other. The surface ratio will be the common area within the two cones on the unit sphere in comparison with the whole surface. Let's choose a Cartesian coordinate system, direct the 'z' axis to the axis of one of the cones, place 'x' axis in the plane determined by the axes of each cone and 'y' perpendicular to them (see FIG. S2). The surface to be calculated Figure S2: Calculation of 1/4 of the surface integral belonging to p eI . The surface to be calculated is limited by thick red lines, the blue colored arrows are the axis of the cones, whereas each plane is a symmetry plane.
posses two planes of symmetry: one is the 'xz' plane at y=0 and the other is the bisecting plane inclined by α I /2 degrees to the z-axis (tan α I /2 = x/z on the plane). Taking that small surface element on the unit sphere (S I ), which is limited by the two symmetry planes and the intersecting cone, the corresponding probability: (S7) In order to calculate S I , we parametrize its surface to spherical coordinates. Applying the Cartesian-spherical conversion form on unit of spheres (x, y, z) = (sin θ cos φ, sin θ sin φ, cos θ), the limits are: Thus, the individual ratio is obtained by evaluating the following surface integral: Introducing the A = cos α I and G = cos γ I notations, the individual surface ratio of the edge becomes Evaluating it to tetrahedral (A = −1/3, G = 0), octahedral (A = 0, G = 1/3) and icosahedral (A = 1/ √ 5, G = 2/3) symmetries, we get: The individual surface ratios should be summed up in general, or simply multiplied with the corresponding number of edges in regular (aka Platonic solids) case to obtain p e , then applying EQ. S5 and S6 to get any of the P i random orientation probabilities (TABLE S1). Note the exact agreement with the calculated value of P 1 of Rey [1] for tetrahedral molecules.

Details of the performed molecular dynamics simulations
For sulfur hexafluoride, calculations at various thermodynamic conditions done, where published diffraction datasets [16,22,23,51] were available. Thus, both NpT and NVT simulations have been performed in the supercritical state (temperature of 398 K and pressures of 128 bar, 155 bar, 394 bar and 1827 bar [23]), in the orientationally disordered/plastic crystalline phase I (190 K and 1 bar [16,51]), in the liquid state close to the triple point (225 K and 10.0552 bar -no measured diffraction data are available here). Only NVT calculations have been performed in the gaseous phase (293.15 K and 10.2 bar [22]). Initial configurations containing 5000 molecules are created for gaseous, liquid and fluid phases by randomly placing and orienting molecules into the simulation box determined by the density. In the crystalline bcc phase, each of the 5488 molecules (cell size 14x14x14 of the crystallographic unit cell with 5.89 Å lattice constant) are oriented randomly around equilibrium positions of sulfur.
For fullerene, only its plastic crystalline phase is modelled at room temperature by NVT simulations using three forcefields (FIG. S2). Each simulation started from an 8x8x8 fcc unit cell (corresponds to 2048 molecules, with lattice constant of 14.16 Å) points as centres of gravity of C 60 , where molecules randomly oriented. Table S2: Lennard-Jones parameters and partial charges applied for fullerene MD simulations at room temperature, ambient conditions. For flexible molecules (applied with 'Girifalco' [41] and 'Lu' [44] forcefields), the flexible intramolecular potential set of Monticelli et al. [42] has been used for carbon atoms: bondlengths of 0.14 and 0.145 nm with 392459.2 kJ mol * nm 2 , bond angles of 120 and 108 degrees with 527.184 kJ mol * rad 2 , improper dihedrals of 143 degrees with 100 kJ mol * rad 2 force constant.

Forcefield
Carbon During the simulations, the Verlet cut-off scheme has been used with van der Waals cut-off at 2.0 nm (1.0 nm for C 60 ). For electrostatics, the particle-mesh Ewald algorithm has been applied, with the same Coulombic cut-off distance when partial charges were non-neutral. The simulation time step was 1 fs. The configuration energies were minimized first, then equilibrated at different levels: for 0.5 ns Berendsen-thermostat (τ T = 0.2 ps, for C 60 it was 2 ps) applied at the right temperature first, then Nose-Hoover thermostat for 2 ns. For NpT runs, an additional 200 ps by Berendsen barostat (with τ p = 1.0) and at least 1 ns by Parinello-Rahman barostat were applied. Carrying out equilibration, production runs were started for 10 ns, saving coordinates at each 100 ps.
In the case of simulations with the 'Lu' and 'Sprik' forcefields of C 60 , a modelling sequence at 700K, then at 500K and finally at 300K has been conducted, where each modelling started from the final configuration of the previous simulation.
For density determination NpT, while for other cases NVT simulations are performed.
6.3. Total scattering powder diffraction patterns: available data and way of calculation from simulations All measured total scattering powder diffraction pattern have been taken from literature [16,22,23,31,51]. Numerical experimental datasets were available for the gaseous [22] and plastic crystalline phase of SF 6 at 190K [16,51] while for supercritical fluid states [23] and C 60 at room temperature [31] datasets have been digitized from figures of differential cross-sections and intensity, respectively.
For crystalline phases, total scattering structure factors have been calculated using a new implementation of the RMCPOW method [52]. For SF 6 11 configurations were used for calculation, separated by at least 1 ns of NVT runs, while only the final configuration was taken for C 60 s. Bragg contributions convoluted by the instrumental resolution function, which was "TOF profile function-2" of GSAS [53] for SF 6 and Gaussian shape for C 60 . Their parameters are optimized by a trial configuration. Finally, individually calculated patterns are averaged.
For other phases, the 'rdf' program of GROMACS provided the partial RDFs, averaged over 101 configurations taken at each 100 ps. Then RDFs are Fourier-transformed and weighted sum created according to neutron coherent scattering lengths and composition.
Generally, during the comparison of the measured and simulated total structure factors, the former is rescaled and a constant background is subtracted to avoid errors from different normalization. For crystalline datasets, the measured pattern rescaled and corrected by 2nd order polynomials above 2.5 Å −1 first, then it has been compared directly over the full measured range.
In order to compare the calculated pattern (I calc ) with the corrected measured one (I meas ), a reliability factor analogous to R wp at Rietveld-refinement [54] is used: The difference is the weights by points are equal and the denominator contains the calculated data instead of the measured one.
6.4. Checking the performance of SF 6 forcefields against total scattering powder diffraction data and density The agreements between measured [22,23,51] and NVT MD simulated structure factors are compared in TABLE S3 using the reliability factor defined in EQ. S15.
Comparing the performance of different forcefields, the worst agreements are found for the 'Strauss' one for each condition, while other forcefields show similar agreement. It is interesting to note that the non-optimized 'Pawley' forcefield shows the best agreement for some states. As a consequence, the 'Dellis', 'Pawley' and 'Strauss' forcefields are selected to represent the corresponding agreement in FIGURE S3 for the disordered phases.
In the gaseous phase (FIGURE S3a), each forcefield provides almost perfect agreements with the measured dataset [22] because the intramolecular structure factor dominates the total structure factor. Focusing on more dense states, the 'Strauss' forcefield total structure factors have not reproduced the measured ones, and close to the triple point liquid phase (FIGURE S3f) large oscillations appeared as the simulated system becomes crystalline. An earlier MD investigation [24] also showed that this forcefield results in peak shifts of the total radial distribution function in comparison with the experimental one, while other forcefields show similar behaviours to each other.
Observing the performance of the other forcefields, agreements at supercritical states (FIGURES S3b-e) are excellent. Only slight differences appear at the low-Q part: the total structure factors of 'Pawley' are slightly less intense than those of 'Dellis': in FIGURE S3d intensity starts to rise below 0.5 Å −1 . The other differences between simulated and measured total structure factors can be explained by errors of digitalization. For the plastic crystalline state (FIGURE S4), the simulated pattern provides satisfactory agreement with the total structure factor. Three main differences can be identified: the Bragg-peak intensities; diffuse scattering intensity differences between 1.5 and 5 Å −1 ; and intensity differences beyond 15 Å −1 . The first one might result from MD simulations provided a slightly disordered model in comparison with reality as the simulated intensity of Bragg-peaks become smaller as Q increases. The second one might come that in this range the measured 3 bank datasets overlap, making their alignment slightly difficult. The last one might originate from the limits of harmonic approximation of intramolecular bond lengths and angle bending instead of a more realistic approximation for far from its equilibrium value. Here, the 'Dellis' model performs better in comparison with the 'Pawley'.  [22,23] total scattering structure factors of gaseous [22] (a), supercritical fluid [23] (b)-(e), and prediction for close to triple point liquid state (f) of SF 6 . Experimental datasets [22,23]: red circles; MD with 'Dellis' forcefield [24]: straight blue lines; with 'Pawley' forcefield [13]: straight orange lines; with 'Strauss' forcefield [23]: straight black lines.
The average densities produced by NpT runs are shown in TABLE S4 for almost every forcefield. It does not contain gaseous phase datasets, because large volume fluctuations prevented performing NpT simulations, while the 'Strauss' forcefield provided a poor agreement to the total structure factors in the NVT ensemble, so the corresponding NpT simulations have not been done. Comparing the experimental with the simulated densities at 6 states, we can conclude that generally, the optimized and the 7 sites 'Dellis' forcefields perform better even for the 190K and 225K states that are below the applied range used for optimization by Dellis and Samios [24]. The strong negative difference in the density of the 'Pawley' forcefield for states near the triple point refers to shifts towards the gaseous state (the difference is notable in terms of total scattering structure factor in FIG. S3f), while positive differences for 'Kinney' refer to some transition to the crystalline state.
Taking into account the agreement with experimental density and with neutron total scattering diffraction patterns, the 7 sites 'Dellis' (with partial charges) or the 6 sites 'Pawley-opt' forcefields are recommended to use. (c) Figure S4: Comparison of MD simulated (with the 'Dellis' forcefield [24]) and experimental [51] total scattering powder structure factors of SF 6 at 190 K in the plastic crystal state. The low-Q-part is shown on (a), while insets show the agreement at (110) Bragg-reflection (b) and high Q (c). Experimental [51] datasets: red circles; MD simulated total (Bragg+diffuse) and only diffuse scattering contributions: straight black and green lines, respectively.