Dynamic simulation of orientational disorder in organic crystals: methyl groups, tri­fluoro­methyl groups and whole mol­ecules

Large amplitude librations of atomic groups or of entire mol­ecules in their crystals are simulated using optimized inter­molecular potentials and crystal structures deposited in the Cambridge Structural Database. The analysis proceeds by a simple static model in which reorientations take place in a fixed environment, or by Monte Carlo (MC) simulation of equilibria dotted by rotational defects, or eventually by full Mol­ecular Dynamics (MD). The simplest approach provides a valuable qualitative preview, but MC and MD are becoming easily accessible to the general solid-state chemist thanks to the facilities of the newly developed Milano Chemistry Mol­ecular Simulation (MiCMoS) platform. Their combined results offer a wealth of information on the behaviour of phen­yl–methyl and phen­yl–tri­fluoro­methyl groups, almost invariably affected by rotational flipping, whose nature and consequences are discussed with respect to disorder modelling in the refinement of X-ray structures. Whole-body reorientation takes place in flat mol­ecules, benzene being the well-known prototype, but also in a very large mol­ecule like coronene. Mol­ecular dynamics of rotations in the cyclo­hexa-1,4-diene crystal offer a spectacular picture of the energetic profiles with jumping times. The dynamic oscillations described here are seldom considered in the formulation of crystal `bonds' or of `synthon' stability.

1. Introduction

Throughout history the notion of `crystal' has been associated with minerals, solids whose constituent chemical units are mostly ionic. When X-ray diffraction came on the scene, the Braggs were able to investigate the inner texture of crystals at the atomic level, and it was soon realized – not without surprise – that most solids are indeed crystalline, that is, endowed with inter­nal long-range periodic symmetry (note how the use of the word `order' has been avoided here). Even more surprising must have been the realization that organic com­pounds are also crystalline, in spite of their being soft and having very low melting points com­pared with salts. The determination of the crystal structures of carbon com­pounds brought, as a revolutionary bonus, the assessment of mol­ecular structures, and for half a century organic X-ray crystallography has been the method of choice for the undisputable determination of mol­ecular shape and size in terms of conformations, chirality and accurate bond lengths.

Intra­molecular curiosity having been amply (one could say overwhelmingly) satisfied, more recently inter­est in organic solids has shifted from mol­ecular properties to materials properties, the ensemble having replaced the constituent as the focus of the scene. In this arena, one cast of mind treats an organic crystal as if it were an organic mol­ecule, stabilized by local inter­molecular bonds between atoms. This simple but popular approach meets an intrinsic difficulty as soon as qu­anti­tative supersedes qualitative – the whole lattice energy of benzoic acid is about one hundredth of the sum of the bond energies within the mol­ecule. Otherwise, for a sound theoretical approach to the organic solid state, one has to stick with fundamental physics. Electronic properties require quantum chemistry, but as concerns thermophysical and structural properties, one is dealing with lattice energies and their derivatives, the lattice forces and, in general, with periodic properties, such as librational and vibrational bands in lattice and mol­ecular dynamics. The first requirement is the preparation of a reliable and manageable model for the inter­molecular inter­action potential among closed shells. Analyses must be based on, and conclusions must be drawn from, the results of com­puter modelling and simulations conducted by means of that numerical puppeteer.

The first part of this contribution is a brief illustration of the derivation and properties of a recent inter­molecular potential energy scheme, preliminary to a second part reporting the analysis of the large amplitude rotational motion (reorientation) of groups or whole mol­ecules within organic crystals at equilibrium. Authors discussing static packing diagrams should consider that at any tem­per­ature, qu­anti­stically even at 0 K, the balls and sticks of those apparently neat drawings are oscillating at a rate of 10¹⁰ to 10¹² times per second. Furthermore, at room tem­per­ature, most organic materials are no more than 50–100 K below their melting points, so that small and com­pact moieties like CH₃ or CF₃ groups, or even entire large flat mol­ecules, may be hopping over some barrier to a com­plete reorientation. Rotational barriers and other dynamic features related to thermal motion in crystals have been reviewed (Dunitz et al., 1988; Trueblood & Dunitz, 1983; Wilson, 2009). The closely related topic of solid-state NMR has also been critically reviewed (Hansen et al., 2013).

Results obtained by simple static models or by Monte Carlo (MC) and Mol­ecular Dynamics (MD) methods are presented. A combination of these methods reveals the detail of reorientation processes and allows a reliable estimation of their timescales, at a fraction of the time and effort needed for an experiment. They also allow some discussion of the more or less obvious consequences that large amplitude motion may have on procedures for the refinement of an X-ray crystal structure determination, especially when the only parameter considered is the R factor. What is unimportant for lowering the R factor may be crucial to the properties of the studied materials.

2. Methodology

2.1. The Lennard–Jones + Coulomb (LJC) potential formulation

The mathematical formulation of an inter­molecular potential function must strike a careful balance between cost and performance. For calculations on single points in phase space-like static lattice energies one can be generous in the physical detail of the potential (Gavezzotti, 2008) or may even adopt an ab initio quantum chemical treatment. But for Monte Carlo (MC) and Mol­ecular Dynamics (MD) simulations, where energies must be calculated 10⁷ to 10¹⁰ times, there is hardly any choice but a simple atom–atom formulation. The LJC potential between atoms in different mol­ecules at a distance R_ij and with point charges q_i reads:

where A₁₂ and A₆ are empirical positive coefficients. The first term is supposed to represent the repulsion arising from Pauli-overlap denial among closed shells and the second term should account for the qu­anti­stic hyperpolarization called `dispersion'. Fig. 1 shows the shape of the LJ 6–12 part. Taken at face value, these curves show that equilibrium in condensed matter comes from a com­promise between dispersion stabilization and destabilizing repulsion. The adherence to physical principles is, however, limited, considering that a simple inverse-distance law must account for such com­plex phe­nom­ena, and that no explicit primary polarization term is present; these closed-shell approximations break down, both physically and mathematically, at very short inter­atomic distances. In recent work at the Milano crystal chemistry group, universal coefficients of the expansion have been calibrated along with point charges of high quality (the MI-LJC scheme) for use with organic crystals and liquids, with excellent performance against a variety of experimental checkpoints (Gavezzotti et al., 2020). The need for this longish exposure, pointing out what seem rather obvious facts, will be apparent in the discussion of the nature of reorientation barriers.

Figure 1 The inter­molecular Lennard–Jones potential exemplified by the Cl⋯Cl potential in the MI-LJC scheme; note the unrealistic divergence of the two separate branches at short distance.

2.2. Data selection: overall survey of crystal disorder

The 2018 Cambridge Structural Database (CSD; Groom & Allen, 2014; Groom et al., 2016) was searched for the text `disorder' with N_at ≤ 35, one residue and Z′ = 1. These restrictions make the number of hits manageable and place a better focus on disorder without causing undue bias or loss of generality. The search yielded 2580 hits: static substituent disorder on polysubstituted aromatic rings, disorder of ali­phatic side chains, head-to-tail disorder in elongated mol­ecules (e.g. azulene), the very frequent oscillation of H atoms in methyl groups and of F atoms in tri­fluoro­methyl groups, and proton hopping in the cyclic hydrogen bonds of carb­oxy­lic acids. Statistics by machine analysis of these occurrences is impossible, because disorder is described with widely different locutions due to the freedom allowed by the CIF format. Coordinates of disordered atoms may be deposited in full, or just for one of the positions used in the refinement, or may be altogether absent for H atoms; in the last two cases, some essential structural information is missing.

Our inter­est here is restricted to large amplitude mol­ecular motions, namely, (a) C—CH₃ groups in which the H atoms rotate around the C—C axis; (b) C—CF₃ groups, as above; (c) flat rigid com­pounds that perform in-plane motions of the entire mol­ecule. The purpose is the estimation of energy barriers to these motions, the study of the distribution of rotation angles at equilibrium and an estimation of dynamic reorientation rates.

2.3. Crystal and mol­ecular data preparation

Structure information is acquired using the appropriate modules in the MiCMoS environment (https://sites.unimi.it/xtal_chem_group/index.php). Data are retrieved from the CSD (Retcif module) and the centre of the atomic coordinates is reset as close as possible to the unit-cell origin, an important precaution when running lattice energy calculations, to generate com­pact coordination spheres.

Retrieved items are then sieved by downstream modules, editing out incom­plete and erroneous entries that have nevertheless passed all the deposition checks. C—H, O—H and N—H distances are renormalized as usual by the Retcif/Retcor module sequence, but the many entries that lack coordinates of hydrogen-bonding H atoms are useless; a com­prehensive study (Gavezzotti, 2021) finds that only about 15% of the CSD entries are suitable for inter­molecular energy analysis of organic com­pounds. Only one set of atomic coordinates is retained for the disordered groups – which set is chosen is immaterial since the energy calculations will anyway sweep the entire set of possible rotational positions. An MP2/6-31G** quantum chemical calculation is carried out to provide the `ESP' charge parameters needed in the MI-LJC scheme. The final result of these preliminaries is a MiCMoS oeh-type file with unit-cell dimensions, atomic coordinates, atomic point charges and symmetry operations in matrix form (a tem­plate is provided in §S-1 in the supporting information). This file type is ready for inter­molecular energy calculations, as well as for starting Monte Carlo and Mol­ecular Dynamics simulations.

2.4. Inter­molecular barriers: zero-level approach

The total barrier for reorientation in the crystal is the sum of an intra­molecular and an inter­molecular contribution. Intra­molecular torsion profiles are obtained by MP2/6-31G** mol­ecular orbital calculations on sample com­pounds, com­puting total energies as a function of the rotation angle: they are 0.1 kJ mol⁻¹ in Ph–CH₃ com­pounds and 0.6 kJ mol⁻¹ in Ph–CF₃, both values being in practice indistinguishable from zero. Ordering in these groups depends exclusively on the inter­molecular field.

The zero-level (ZL) approach to the simulation of inter­molecular energies for rotational mol­ecular rearrangements in crystals is as follows. A crystal slab is built by ±2 cell translations in the three directions, yielding a central reference mol­ecule (RM) plus 125 Z − 1 surrounding mol­ecules (SM) (Z is the number of mol­ecules in the unit cell). The coordinates of some or all of the atoms in the RM are rotated around the chosen axis; the matrix algebra (Gavezzotti & Simonetta, 1975) is specified in the supporting information (§S-2) and the program code with input–output examples is deposited as supporting information (§S-3). The total potential energy of the system is calculated by summing all atom–atom contributions of Equation (1) between the RM and all the motionless SMs in the crystal slab as a function of the local rotation angle, yielding an energy profile for the mol­ecular rearrangement in a fixed crystal environment. This is the energy needed to introduce one mole of single orientational defects within the crystal. In this approach, the energy barriers are presumably an upper-limit estimate, because intuitively, correlation with a flexible environment is more likely to offer com­pliance than resistance to the com­pressive rotation. (Note that the 1975 article used lattice energies instead of potential energies, thus yielding halved barrier values.)

2.5. Monte Carlo approach

A special feature of the Monte Carlo (MC) Mcmain module of the MiCMoS environment, specifically oriented to small-mol­ecule condensed phases, is that some parts of the mol­ecule may be kept rigid while some parts are allowed torsional freedom. In this way, mol­ecules are assigned six overall rigid-body degrees of freedom (dof) plus one dof for each allowed torsion; for example, tri­fluoro­methyl­benzene is a 7-dof system. Thus, no simulation time is spent in sampling irrelevant parts of intra­molecular phase space. A com­putational box for the crystal model is set up by an appropriate number of repetitions of the unit cell along the three directions of space (MiCMoS module Boxcry), and module Pretop provides a tem­plate for the force field `topology' file. In our case, the only allowed intra­molecular degree of freedom is the torsional CH₃ or CF₃ motion, described by the usual (Gavezzotti & Lo Presti, 2019) cosine potential functions (see also the supporting information, passim). Ordinary NPT runs (constant number of particles, pressure and tem­per­ature) with periodic boundary conditions are then carried out under the MI-LJC inter­molecular potential until a steady state (equilibration) is reached, after which the distribution of torsion angles or of whole-mol­ecule orientations is analyzed. MC does not provide energy profiles or barrier heights, but its advantage over ZL is that its final frame reproduces the equilibrium distribution of rotational defects. More detail can be found in the extensive documentation and tutorials on the freely available MiCMoS site.

2.6. Mol­ecular dynamics approach

The preliminaries to the mol­ecular dynamics (MD) approach are identical to those for the MC approach, but MD samples the entire phase space, with potentials and forces com­puted over all stretching, bending, torsion and inter­molecular degrees of freedom. After the Boxcry and Pretop stages, unconstrained dynamic simulation in all the 3N − 6 degrees of freedom is carried out in module Mdmain, with periodic boundary conditions, tem­per­ature and anisotropic pressure control, for the time needed in order to observe possible rotational jumps. Intra­molecular bond stretching, bond bending and torsional degrees of freedom are taken care of by the standard MiCMoS force field embedded in Pretop (Gavezzotti & Lo Presti, 2019; see also the supporting Information, passim). Dynamic events are monitored by following in time the distribution of relevant torsion or rotation angles, so that trajectory analysis provides a mol­ecular level picture of the proceedings and an estimate of the time lapse during and between rotational jumps.

2.7. Reproducibility

Details of all mol­ecular models, force fields and MC and MD operating conditions with input–output examples are deposited in the supporting information (§S-4). Further details can be found in the documentation on the MiCMoS site. Together with the public availability of the MiCMoS codes, this ensures com­plete reproducibility of all calculations carried out in this work. Using a workstation or even a moderately powerful laptop com­puter, typical com­puting times for the described simulations are a few seconds for ZL, one hour for MC and a few hours for MD.

3. Results and discussion

3.1. Zero-level simulations: methyl-group rotation in methyl­benzenes

A sample of 29 crystals of com­pounds with ordered methyl groups attached to a phenyl ring without ortho substituents or other intra­molecular crowding conditions was selected; five more structures with methyl groups declared as disordered were added (see ref­codes in §S-5 of the supporting information). Aliphatic methyl groups were not considered because the high intra­molecular barrier (≃15 kJ mol⁻¹) prevents the rotation. In the ZL approach, the C(Ph)—C(Me) bond is set as the rotation axis for the three H atoms, and the rotational energy profile is scanned in steps of 5° and is characterized by: (1) E°, the relative energy at the position of the starting coordinates; (2) θ_min, the C=C—C—H rotation angle at which the energy is minimum; (3) K, the barrier, that is, the highest relative energy at rotation angle θ_max.

Fig. 2 shows the extremes: in the ideal situation, E° and θ_min are both close to zero, the barrier is relatively high and perfectly threefold periodic with θ_max ≃ 60 ± 120°; for a disordered case, the H atoms have been placed in a position that does not correspond to a minimum energy, although differences of less than 1 kJ mol⁻¹ are scarcely significant. Fig. 3 shows a landscape of methyl rotation barriers, whose top is generally near 60° rotation, as it should, halfway between two minima, with very nearly zero relative energy (E°) at the experimental position (θ_min = 0°). A few points at E° = 2–4 kJ mol⁻¹ and θ_max = 0–20° suggest that some H atoms (mostly those introduced to model the disorder) have been placed at inaccurate positions in the experimental data. The barriers to rotation are very low in four of the five structures labelled as disordered, but otherwise the many calculated barriers of 3 kJ mol⁻¹ or less suggest that disorder may have been frequently overlooked. As might have been guessed, the origin of the rotation barrier is invariably atom–atom repulsion, with insignificant variation in Coulombic or dispersion energies. The barrier height depends on the details of the inter­molecular environment; for example, the structural reason for the very high barrier in one outlier in Fig. 3 can be traced to an inter­locking of nearest-neighbour methyl H atoms opposing H-atom mobility.

Figure 2 The ideal energy profile (blue line) for rotation around the C—Me axis of a Ph–Me methyl group in an ordered structure (CSD ref­code BAZFUG), and the wiggly profile for a disordered case (CSD ref­code GOFNEW). These are total barriers since the intra­molecular contribution is virtually zero.

Figure 3 The zero-rotation energy, E° (top), and the barrier to rotation for 29 ordered and five disordered Ph–CH3 groups (bottom). θmax is the highest-barrier rotation away from the X-ray position (ideal value 60°), the corresponding θmin is in all cases equal to θmax ± 60°. Filled red circles are structures reported with methyl-group disorder. The inset shows the inter­locking of methyl H atoms (magenta) that produces a high barrier in the outlier (CSD ref­code OMONEJ). The axial label is the same in both frames.

3.2. Zero-level simulation of CF₃ rotation in tri­fluoro­methyl­benzenes

The analysis follows the same lines as for methyl com­pounds. Rotation energy profiles are more wiggly, barriers are much higher than for methyl-group rotation and 120° periodicity is not always strictly observed. These are effects of the size difference of the rotating group: a bulkier and more extended object (C—H = 1.08 Å and C—F = 1.35 Å; H-atom radius = 1.10 Å and F-atom radius = 1.45 Å) has more chances to `bump' into its surroundings. Fig. 4 shows the distribution of rotation barriers in Ph–CF₃ groups (ref­codes are available in §S-6 of the dupporting information): the trend is rather neat, with barriers for disordered groups nearly all below the 20 kJ mol⁻¹ line. A further confirmation of the reliability of the indication from ZL calculations comes from the cases shown in Fig. 5, where the same mol­ecule has one ordered and one disordered CF₃ group, with a widely different rotation barrier due to a clear difference in the inter­molecular environment.

Figure 4 The zero-rotation energy, E° (top), and the barriers to rotation for a sample of 35 disordered (full red circles) and 20 ordered (black squares) Ph–CF3 groups (bottom). θmax is the highest-barrier rotation away from the X-ray position (ideal value 60°). The axial label is the same in both frames.

Figure 5 Compounds with one CF3 group free (disordered) and one CF3 group inter­locked with the surroundings (ordered): CSD ref­codes ENEHUB (top) and EDOROG01 (bottom). The red dashes denote the inter­locking contact. The numbers in parentheses are the calculated ZL barriers (in kJ mol^−1).

3.3. ZL simulation of whole-body rotations

Mol­ecules whose external envelope is disk-like or globular may sometimes undergo whole-body reorientations in their crystals. A typical example is benzene (ref­code series BENZENXX; see appendix A), along with its sixfold expansion coronene (ref­code CORONE02). Fig. 6 shows the results of the ZL com­putational experiments. The barrier is surprisingly low, even for coronene, and reorientation is the reason for the unusually high melting points of these two com­pounds (706 K for coronene) due to low melting entropy. Fig. 6 also shows the energy com­ponents of the com­putational barrier for coronene: Coulombic terms are insignificant, in spite of their often invoked role in so-called `C—H⋯π' inter­actions. At the top of the barrier, some atoms are in short inter­molecular contact with static neighbouring mol­ecules; this would induce a dispersive stabilization, but repulsion rises faster and the net result is destabilization. This explanation stems from the shape of the curves in Fig. 1. Empirical curves reproduce physical effects, at least when not far from their minima, so this result demonstrates how mol­ecules in crystals are positioned at a point of balance between the requests of attractive and repulsive forces.

Figure 6 (Top) The ZL energy profiles for in-plane rotation of benzene (red) and coronene (black, inset) in their crystal. The secondary minimum at 30° for coronene corresponds to exposure of the hydrogen-free bay area. (Bottom) The com­ponents of the barrier for coronene: stabilizing dispersion (blue dots) against increasing repulsion (black broken line). The thick red line denotes insignificant Coulombic terms. The 120–360° landscape is periodic.

3.4. Temperature effects

The restriction to group rotation has an intrinsic structural explanation due to inter­molecular blocking, but also density variations with tem­per­ature should have an influence. Fig. 7 shows the distribution of ZL rotation barriers against the tem­per­ature of the X-ray work. There is no discernible trend, but at room tem­per­ature, all methyl rotation barriers are below 5 kJ mol⁻¹, while most of the higher CF₃ barriers for ordered structures are at low tem­per­ature. In order to observe static ordered Ph–CH₃ or Ph–CF₃ groups, X-ray determinations should be carried out at low tem­per­ature and, vice versa, any crystal containing these groups is likely to have them disordered not too far from the melting point.

Figure 7 Rotation barriers against tem­per­ature of the X-ray determination for Ph–CH3 and Ph–CF3 groups.

Fig. 8 shows the effect of tem­per­ature and pressure on the height of the rotational barrier in benzene, as obtained by ZL calculations on crystal structures determined at various tem­per­atures and at high pressure. Both trends are rationalized on the basis of variations in density and in inter­molecular com­pression. The experimental tem­per­ature of the rotation onset and the rotation barrier, as determined by solid-state NMR, are around 120 K and 15–17 kJ mol⁻¹, respectively (Wendt & Noack, 1974; Andrew & Eades, 1953). These two experimental data are in almost qu­anti­tative agreement with the com­putational results shown in Fig. 8.

Figure 8 The ZL com­putational barrier to rotation of benzene in crystal structures determined at various tem­per­atures. The red square at the top right is for the structure under pressure.

Extensive calculations carried out on crystals of many flat benzene derivatives show that attached atoms other than hydrogen prevent rotational jumps. For example, none of the fluoro- or chloro­benzene crystals, however substituted, show a ZL rotational barrier below 50 kJ mol⁻¹. For hexa­chloro­benzene, the barrier is 96 kJ mol⁻¹ at 100 K, and is down to 61 kJ mol⁻¹ at room tem­per­ature. This crystal probably becomes rotationally disordered at higher tem­per­ature, its melting point being also unusually high, at 504 K. Another can­didate for in-plane rotation is cyclo­hexa-1,4-diene, with a ZL barrier of only 25 kJ mol⁻¹ (see below).

3.5. Monte Carlo simulations

Standard MC simulations were carried out for some of the crystals previously studied by ZL calculations (more detail is given in §S-7 of the supporting information). MC provides a picture of the equilibrium inter­nal structure of a bulk phase, including rotational defects; for the present purpose, the relevant parameter is the distribution of torsion angles of selected groups or the spread of mol­ecular orientations. Fig. 9 shows some results: a structure with ordered CF₃ groups and a very high ZL barrier (67 kJ mol⁻¹) has a narrow peak around the experimental torsion angle, while a structure with disorder and a very low barrier (6 kJ mol⁻¹) shows a spread over the whole range. Notably, even in the ordered structure, the distribution peak has a width at half maximum of about 50°, hinting at significant oscillation freedom. A tentative explanation of the structural reason for the barrier difference in terms of inter­molecular environment is given in §S-8 of the supporting information.

Figure 9 Distribution of C—C—C—F torsion angles from a Monte Carlo simulation of a bulk crystal phase with ordered (black line; CSD ref­code EFITAO) and disordered (red squares; CSD ref­code VILBUN) CF3 groups at room tem­per­ature. Zero is the experimental position for the ordered phase, or for one of the two sets of F-atom positions in the disordered crystal.

In a more revealing com­puter experiment, the crystal structure of N,N-dimethyl-N′-[3-(tri­fluoro­meth­yl)phen­yl]urea, reported with disordered CF₃ groups at room tem­per­ature (Yu et al., 2008), was also simulated at the com­putational tem­per­ature of 100 K. The result is summarized in Fig. 10; at room tem­per­ature there is a significant distribution of torsions between −30 and +30°, while at 100 K the rotational freedom is quenched, with zero population in the ±15° range. The MC simulation thus predicts an ordered crystal structure at low tem­per­ature.

Figure 10 The distribution of C—C—C—F torsion angles in a crystal structure disordered at room tem­per­ature (CSD ref­code HODHIS) and in the same structure at 100 K simulated by Monte Carlo. Mol­ecules related by centrosymmetry have torsions of opposite signs, hence the doublet in the profile.

For mol­ecules of appropriate shape, inter­nal rotation in organic crystals borders with the science of liquid or plastic crystals for materials like the adamantanes (Brand et al., 2002) or flat circular polyacenes (Zhong et al., 2018). Coronene is an example of the latter category, having a ZL rotation barrier of only 26 kJ mol⁻¹ (Fig. 6). Monte Carlo simulations of the coronene crystal have been carried out at 300, 400, 500 and 700 K, monitoring the distribution of distances between centres of mass with their corresponding angles between in-plane vectors. The results (Fig. 11) show that some reorientation occurs already at room tem­per­ature; the orientation spread becomes diffuse at 500 K, while at the melting tem­per­ature (see §S-9 in the supporting information), there is also an incipient randomization of the distances between centres of mass, which is proper for a crystal that is about to lose both long- and short-range periodicity on the way to the melting transition. Incidentally, the MC simulation provides an estimate of the thermal expansion coefficient of the crystal, (1/V) dV/dT = 7 × 10⁻⁵ K⁻¹, which is rather low but in line with normal values for organic crystals.

Figure 11 MC simulation of the coronene crystal (P21/n, Z = 2): distribution of mutual in-plane rotation angles against distances between mol­ecular centres of mass (ordered crystal values: 0° at 4.69 and 9.39 Å; 132° at 8.42 Å). Restricted distribution with some 60° jumps at room tem­per­ature (blue diamonds) and spread in distance and angle at 500 K (red circles).

3.6. Mol­ecular dynamics simulations

An obvious candidate to probe the dynamics of mol­ecular librations in tri­fluoro­methyl groups is the parent com­pound tri­fluoro­methyl­benzene (TFMB; CSD ref­code XOGJAG). Its ordered crystal structure has been determined at 213 K (Merz et al., 2014), just below its melting tem­per­ature of 242 K. The calculated ZL barrier to CF₃ rotation is 26 kJ mol⁻¹, in the correct zone for the population shown in Fig. 4. For com­parison, two other structures, both determined at room tem­per­ature, were considered, i.e. 5-(3-tri­fluoro­methyl­ben­zyl­idene)thia­zolidine-2,4-dione (Bruno et al., 2002) (TFTD; m.p. 453 K; CSD ref­code EFITAO; see scheme in Fig. 9, top left), with a very high ZL rotational barrier of 67 kJ mol⁻¹, and tri­fluoro­methyl­phenyl­propanoic acid (Guan et al., 2009) (TPPA; m.p. 379 K; CSD ref­code VOQLUJ), with a very small ZL barrier of 6 kJ mol⁻¹.

A mol­ecular dynamics simulation was carried out for each of the three com­pounds, with standard conditions (see §S-10 of the supporting information for more detail), applying the C—C—C—F intra­molecular torsional potential and keeping the CF₃ group rigid by the usual MiCMoS device of stiff C—F stretching and C—F—C bending potentials. Fig. 12 shows that in a very short time lapse, the distribution of the torsion angles for TPPA randomizes as expected, indicating widespread reorientation, and that, rather surprisingly, the same happens to a smaller extent for TFMB. The same plot for TFTD gives no spread even after a 20 ps MD run, confirming the MC result of Fig. 9.

Figure 12 Distribution of torsion angles after a 10 ps MD simulation of the crystal structure of (top) TPPA at room tem­per­ature and (below) TFMB at 213 K. The angles are ±30 and ±50°, respectively, in the static crystal structures (zero time). The analogous plot for TFTD shows no spread after 20 ps. Axial labels are the same in both frames.

The time monitoring of torsional libration is performed by preparing a trajectory of rotation angles using successive MD frames at steps of 0.1 ps (see §S-11 of the supporting information for details of this construction). The trajectory is then scanned to find the minimum and maximum oscillation values, and their difference (the oscillation amplitude) is plotted against the time lapse between the two extremes. The result is shown in Fig. 13. For TFTD, far from the melting point and with a high ZL barrier, the average oscillation is restricted and sluggish, becoming smaller at higher time spans, quite in keeping with the ordered structure found in the X-ray analysis. For TFMB, close to the melting point, the average oscillation range about the equilibrium position for most mol­ecules over a time span of a few picoseconds is wider, but about 25% of the mol­ecules perform an oscillation larger than 120°, thus apparently approaching or even overshooting the top of the barrier at ±60°. The part of the plot related to TPPA shows com­plete rotational freedom as CF₃ groups perform entire or even double barrier jumps in a very short time span.

Figure 13 Range of torsional oscillation of CF3 groups in an MD simulation of the TFMB (213 K), TFTD and TPPA (295 K) crystals as a function of the time lapse between the two extremes. Trajectories are sampled at 0.1 ps inter­vals.

No disorder was reported in the crystal structure of TFMB although the anisotropic displacement parameters (ADPs) of F atoms show a circular anisotropy hinting at rotational freedom (the R factor is 7.99%). In view of the results shown in Figs. 12 and 13, however, the discussion of fine features of C—H⋯F or F⋯F contacts and relative stabilizations seems at least controversial, with atom positions oscillating by ≃1 Å (≃40° on an 8 Å circumference). Neglect of thermal motion is one of the main (if seldom pointed out) shortcomings of all approaches based on localized atom–atom bonds with related crystal engineering exertions. The structure of TPPA (R factor 6.8%) was modelled by a six-site CF₃ group; the total rotational freedom could probably be better modelled by a continuous ring of F-atom density, something that cannot be done with standard X-ray data treatment packages.

As guessed earlier, cyclo­hexa-1,4-diene (Jeffrey et al., 1988; CSD ref­code VACCEH, 153 K, melting tem­per­ature 223 K) is a candidate for in-plane reorientation. A 100 ps MD run was carried out for this crystal, and the distance–orientation analysis on the final frame confirms com­plete rotational disordering at its melting tem­per­ature (see §S-12 in the supporting information). The time analysis of the dynamic treatment offers a firsthand view of the mol­ecular proceedings. Fig. 14 shows two typical paths traced by a sample mol­ecule: the first has three sudden 60° jumps to com­plete inversion of orientation and the second shows a sudden 60° jump with immediate back rotation, before an 80° jump that, being incom­patible with the approximate sixfold symmetry of the mol­ecule, requires a rise in energy. We believe that such a result is an impressive illustration of the power of MD simulation in tracing not only the structural variation, but also the timescale of the portrayed events, offering a significant input to the inter­pretation of X-ray and solid-state NMR experiments.

Figure 14 Mol­ecular dynamics simulation of the cyclo­hexa-1,4-diene crystal (at 223 K), showing the simultaneous structure–energy timeline for two sample mol­ecules. Zeroes are the X-ray position and its potential energy. In the top frame, the energy toll never exceeds 16 kJ mol^−1, much less than the ZL barrier of 26 kJ mol^−1. The vertical axis is in degrees (°) for angles (circles) and kJ mol^−1 for energies (red diamonds). Axial labels and units are the same in both frames.

4. Summary and conclusion

A general survey of CSD entries labelled as disordered has been conducted. Two sets of crystal structures for phen­yl–methyl and phen­yl–tri­fluoro­methyl com­pounds have been picked from this landscape. In a first approximation, the inter­molecular potential energy profiles for the rotation of CH₃ or CF₃ groups around the Ph—C axis in one central mol­ecule, surrounded by a rigid environment of motionless neighbours, have been evaluated. The zero-level (ZL) barriers thus obtained are total barriers, since the intra­molecular contribution is negligible. They are upper limit estimates, as there is no com­pliance of the medium around the rotational defect, yet they offer a reliable preliminary indication of the amount of restriction to rotational oscillation. We find that at high enough tem­per­ature (typically not far from the melting tem­per­ature, and thus often even under room conditions) the methyl groups are almost invariably rotating freely over barriers of the order of RT, and the CF₃ groups are also very frequently prone to rotational disorder. A case-by-case check of the barrier height against the treatment of disorder in the refinement of the X-ray crystal structure determination is impossible for the large number of structures in our samples; from sample observations one gets the impression that there is very little preoccupation for structure and properties, the only aim being the placement of fractional atoms somewhere to lower the R factor slightly. For example, in disordered CF₃ groups deposited with multiple F-atom locations, one sees C—F distances varying from 1.20 to 1.40 Å, with C—C—F and C—F—F angles often at impossible values; the best description would be a torus of 3-electron or 27-electron density at the proper distance from the carrier C atom. This issue borders on the sociology of X-ray diffraction work: no one is inter­ested in the geometry of peripheral groups if the determination is a routine one. The point we want to make here is that if one is inter­ested in materials properties, the possible disorder can be easily investigated by ZL-type calculations, requiring a few minutes of file editing for input to the program Cryrot (see §S-3 in the supporting information) that runs in a few seconds.

ZL calculations for the reorientation of entire mol­ecules for flat rigid com­pounds have been also carried out. The main result is that the possibility of rotational disorder is very restricted by a shape factor; besides the obvious and well-known case of benzene, for example, there is easy rotation in a large mol­ecule like coronene, but rotation is blocked in a simple mono-derivative like fluoro­benzene. Hexa­chloro­ben­zene is a candidate with rotational disorder near its melting point.

A more proper description of disorder is by Monte Carlo simulation, yielding a reliable picture of the distribution of rotational defects within the crystal at equilibrium, with the added bonus that one is not restricted to the tem­per­ature of the X-ray determination. The com­putational effort for a Monte Carlo simulation is one order of magnitude greater than for a ZL simulation, but the retrieved information is also one order of magnitude more relevant for materials science. In a rather simple com­putational experiment, a tri­fluoro­methyl group disordered in the experimental structure at room tem­per­ature is predicted to be ordered at 100 K. For benzene, there is excellent agreement between calculation and experimental values for the rotation onset tem­per­ature and the barrier height. The coronene crystal structure investigated between 300 and 700 K shows the ever increasing amount of rotational defects already present at room tem­per­ature. Aside from its value for the general progress of organic solid-state chemistry, in principle, this kind of information could be of great practical value in the planning of materials or of devices based on them.

The ultimate tool in the description of mol­ecular events is Mol­ecular Dynamics simulation, including all degrees of freedom and, most important of all, the time coordinate. A com­parative study of CF₃ motion was carried out on a disordered structure (TPPA) and an ordered structure (TFTD). The final frames show the expected spread for the former and restricted rotation for the latter, but the relevant information is that the rotational jumps are very fast, with the disordered groups performing a few full rotations in less than 10 ps. The parent com­pound tri­fluoro­methyl­benzene is a borderline case, with inter­mediate behaviour. The MD simulation reveals rotational motion in the determined crystal structure with oscillation of F-atom positions by as much as 1 Å.

An example of the most com­plete information provided by MD is a study of the in-plane rotation in the crystal of cyclo­hexa-1,4-diene, whose crystal structure was presented in a short note with the aim of elucidating the mol­ecular structure without consideration of inter­molecular or dynamic facts. The com­puter simulation provides a timeline with orientation and relative energy for each mol­ecule, highlighting the timescale and the amplitude of the rotational jumps: a mol­ecule is seen performing three instantaneous 60° rotations, remaining in each of these configurations of 20–40 ps, on its way to a com­plete inversion of orientation. This result would, among other things, be a perfect com­plement to solid-state NMR experiments.

Of course, all the described conclusions depend on the quality of the applied force field and on a number of com­putational approximations. The MI-LJC potential scheme and the MiCMoS platform in which it is embedded have passed a large number of validation tests, including that presented in this article with the reproduction of tem­per­ature onset and barrier height for benzene reorientation.

If asked what is the most important conclusion of this work, the author would say that such results should be kept in mind when discussing fine detail of atom–atom inter­molecular `bonds', or `short' distances, based on the averaged positional picture provided by X-ray diffraction. Many of the postulated crystal bonds and of the related crystal engineering inferences depend on ghost atom locations that come from plain and often uncritically accepted artifacts of X-ray diffraction structure refinements.

APPENDIX A

Literature citations for structure determinations not explicitly mentioned in the text, labelled by the corresponding ref­code: EDOROG01 (Abe et al., 2012); OMONEJ (Gdaniec et al., 2004); VILBUN (Grunewald et al., 1991); ENEHUB (Jonet et al., 2011); GOFNEW (Kumar & Perumal, 2014); BAZFUG (Xiong et al., 2017); CORONE02 (Krygowski et al., 1996). The benzene crystal structures are: BENZEN06 (15 K) and BENZEN07 (123 K) (Jeffrey et al., 1987); BENZEN20 (100 K) (Woińska et al., 2016); BENZEN18 (150 K) (Nayak et al., 2010); BENZEN25 (225 K) and BENZEN26 (270 K) (Bujak & Mitzel, 2018); BENZEN11 (high pressure; Budzianowski & Katrusiak, 2006).