Discovering structure–property relationships for the phonon band structures of hydrocarbon-based organic semiconductor crystals: the instructive case of acenes

By studying the low-frequency phonon bands of a series of crystalline acenes, this article lays the foundation for the development of structure–property relationships for phonons in organic semiconductors. Combining state-of-the art quantum–mechanical simulations with simple classical models, we explain how and why phonon frequencies and group velocities do or do not change when varying the molecular and crystal structures of the materials.


Electronic Supplementary Information
Discovering Structure-to-Property Relations in Phonon Band

Structures of Hydrocarbon-Based Organic Semiconductor Crystals:
The Instructive Case of Acenes Tomas Kamencek and Egbert Zojer*

S1. Data availability
The simulation input and output required to reproduce the presented results are publicly available in the NOMAD database. 1

S2. Experimental starting geometries
The crystal structures used as starting points for the geometry optimisations are listed in Table   S

S3. Details on the structural similarities and differences among the discussed systems
Table S 2 contains the main parameters of the unit cells of the studied systems. One observes rather significant variations of the area of the basal plane, the aspect ratios, and also the nearestneighbour distances as well as the distances to the nearest symmetry-equivalent molecules (which are given by a1 and a2, respectively). Also, some other variations stand out, like the particularly small (large) area of the basal plane for monoclinic benzene (tetracene) and the generally smaller aspect ratios for the triclinic systems (which also result in smaller nearest-neighbour distances), but by and large, the variations in the geometric parameters appear not to be particularly systematic. Thus, in the following table, several geometric parameters are collected that characterise the arrangement of neighbouring molecules in the unit cell. Note that in order to be consistent with the unit cells of the longer acenes (see caption of Table S 2), the long/short molecular axis as well as the lattice vectors must be reassigned for 1A, 4A, and 5A-II. This reassignment is explained for the most complex case of 1A by comparing the crystal structure of 1A and 2A in Figure S 3.  Table S 3  (i) regarding the tilts of the molecules relative to the normal to the a1,a2 plane, there are two groups of acenes: in 1A-3A, the two molecules per unit cell are tilted in different directions, which are symmetrically arranged left above and below the (negative) a1 axis (see 6 th column in Table S 3 and Figure S 5). As a result, the angle between the long axes of the molecules is rather large but continuously decreases from 1A (39.4°) to 3A (13.8°).
Conversely, for 4A, 5A, and 5A-II, both molecules tilt essentially in the same direction, which is nearly diagonal with respect to the a1 and a2 axes. Consequently, the two molecules are no longer symmetry-equivalent (see Figure S 5) consistent with the triclinic structure of the unit cells.
(ii) This different mode of tilting also results in smaller tilt angles in 4A, 5A, and 5A-II compared to 2A and 3A.
(iii) Monoclinic benzene displays a few additional deviations from the other (short acenes): its tilt angle, , is much smaller than that of 2A and 3A, and, more importantly, here the two molecules tilt nearly perpendicular to the a1 axis (i.e., parallel to the a2 axis) (iv) The angle between the short molecular axes in 1A is much larger than in the longer acenes.
Consequently, in 1A the planes of the two molecules per unit cell come to lie close to perpendicular, while they are much less inclined relative to each other in the other acenes.
(v) Finally, two systems also stand out regarding the displacements of equivalent molecules in consecutive layers (see Figure S 4(b)): in benzene, the normal distance between the planes of the molecules amounts to 2.85 Å; i.e., there is a particularly large perpendicular offset, perp, between molecules in consecutive layers, when viewing the structure along the long molecular axis. Conversely, equivalent pentacene molecules in consecutive layers are close to coplanar in 5A-II.

(vi)
The values of perp and short also motivate, why the structure of 5A is more consistent with the shorter acenes than that of 5A-II.
Overall, these data show that there are two qualitatively different groups of systems (1A to 3A vs. 4A, 5A and 5A-II), which within the groups share common geometrical features. Moreover, concerning most quantitative parameters, the structure of benzene quite notably differs from that of the longer acenes in its group.    Figure S 4(a)). The orientation of 1 and 2 relative to the (negative) a1 axis is also summarized in the 6 th column of the table. For benzene, the definition of the long molecular axis is a priori not unique, as illustrated in Figure S 2, but a detailed comparison between the structures of benzene and naphthalene reveals that the green arrow in Figure S 2 denotes the most consistent definition of the long molecular axis (see Figure   S 3) . long and short denote the angles between the long and short molecular axes of the two molecules per unit cell. perp and short are the components of the vectors connecting equivalent molecules in consecutive layers (i) perpendicular to the molecular planes and (ii) within the planes parallel to the short molecular axes, as sketched in Figure S 4 * The sign of the angles relative to -a1 varies here because of the conceptionally different choices of the a3 axes in the CSD for 4A and 5A-II compared to the other crystals, as indicated in Figure  S 1; however, this does not imply a conceptual change of the structure. Here, positive angles refer to anticlockwise rotations.   The discrete k-meshes of electronic wave vectors were adjusted for each system individually to achieve homogenous sampling of the first Brillouin zones, with the density and the k-mesh being chosen such that the total energies were converged below ~0.5 meV/atom. In these convergence tests, an energy cutoff of 900 eV in all systems in combinations with the k-meshes listed in Table   S 4 turned out to achieve the desired level of convergence. These settings were shown to yield high accuracy results in phonon-related properties compared to experimental measurements in crystalline naphthalene 14 .

System (-centred) k-mesh
The headers of the PAW-pseudopotentials that we used for the DFT calculations (for carbon and hydrogen atoms) are listed in Table S 5.

S4.2 Geometry and unit cell optimisations
The atomic positions and the lattice parameters were optimized to residual forces below 1 meV/Å employing the conjugate gradient algorithm within the following procedure: the independent lattice parameters and the atomic positions were optimised at a few fixed unit-cell volumes (using the VASP optimisation keyword ISIF=4). This procedure avoids Pulay stresses 15 to a large extent. 11 The optimum unit-cell volume was subsequently determined by fitting total energies of those structures with optimized lattice constants and atomic positions but fixed unitcell volumes to a Rose-Vinet equation of state 16 . Finally, a last optimization was carried out at the optimal volume extracted from the fitted equation of state to obtain the final atomic positions and lattice parameters.

S4.3 Phonon calculations
Regarding the phonon calculations using the Phonopy code, 17 a displacement distance of 0.01 Å was used to displace the atoms in supercells of system-specifically adjusted extents (see Table S 6) which were chosen based on thorough convergence tests for 2A. 14 For band structure plots, the reciprocal space between each pair of high-symmetry points was sampled by 200 intermediate points, while for quantities that rely on a homogeneous sampling of the entire first Brillouin zone (density of states/group velocities, heat capacity) system-specific meshes of phonon wave vectors, q, were used (see Table S 6). For the group velocities, the meshes were shifted such that they do not include the centre of the first Brillouin zones, . Table S 6 shows two choices of q-meshes which sample the respective first Brillouin zones with different densities. The "dense" meshes (with a maximum integer of 15) were used for the (computationally more demanding) twodimensional densities of frequencies and group velocities, while the "highly dense" meshes (with a maximum integer of 30) were used for densities of states (DOSs), densities of group velocities (DOGVs) and for the calculation of the heat capacity, CV.

S4.4 Densities of states per frequency and per group velocity
Based on the uniformly sampled group velocities and frequencies, a two-dimensional density of states, g, per frequency, , and per group velocities, vg, was calculated using the following definition: Here, 3N is the number of bands, Nq the number of wave vectors used in the sum over the joint index  containing the band index and the wave vector. The delta distributions in Equation (S1) were replaced by Lorentzian functions of finite widths of vg = 0.02 THzÅ (=0.002 kms -1 ) for group velocities and  = 0.05 THz for frequencies (with x representing either group velocities or frequencies): In other words, the finite broadening parameter, x, corresponds to the half width at half maximum of the artificially broadened peaks.
From the density g in frequency-group velocity space, one can, in principle, directly calculate the density of states (DOS; density of states per frequency) and the density of states per group velocity (DOGV) by integrating Equation (S1) over the (continuous) group velocity or frequency variable, respectively: Although these integrations could be carried out numerically to obtain the DOS and the DOGV from g, these two quantities were calculated separately using denser meshes of (phonon) wave vectors, q, using Equation (S3) and (S4) directly. The reason why less dense meshes were used for the computation of the density g in the first place, was a notable decrease in computational efficiency with the density of the q-mesh: for denser meshes, no visually perceivable changes in g were observed any more, while the computational time (due to fact that large two-dimensional arrays had to be stored and changed) significantly increased.
For the sake of completeness, the following plots show the DOS, the DOGV and the twodimensional densities, g, for all studied molecular crystals.  Figure S 8: Two-dimensional densities in group velocity-frequency-space, g (see Equation (S1)), for all studied organic semiconductor crystals. The density is plotted as the (logarithmic) colour scale.

S4.5 Participation ratios
The participation ratio (PR) of a phonon mode  describes the degree of localization of the atomic motion associated with the phonon mode. [18][19][20][21] Therefore, the PR is a convenient measure to quantify the extent of moving atoms in the unit cell. According to its definition, the PR depends on the atomic masses, mi, of the N atoms in the unit cell and on the  th Cartesian component of the phonon (polarisation) eigenvectors, e  ,. To calculate the PR of a mode, one must carry out the above summation over Cartesian directions, , and atoms, i.
One can easily verify that the PR ranges between 1 (if all atoms move with the same amplitude) and 1/N (if only one atom in the unit cell moves).

S5.1 -phonons
To facilitate the association of the provided -phonon animations with the discussed vibrations, the following tables list the frequencies and types of modes that the video files show.

S5.2 Acoustic phonons close to 
In addition to the animations of -phonons, also the acoustic phonons in selected directions were animated to gain some insight into the atomic displacements of the associated motions. The atomic geometries displaced along the phonon eigendisplacements (at q ≠ 0) were generated in a way analogous to the implementation in Phonopy. 17 The so-created displaced geometries were subsequently rendered using the Ovito software 22 and saved as files "TA1.mp4", "TA2.mp4", and "LA.mp4" (for the first transverse, the second transverse, and the longitudinal acoustic modes) in the respective subfolders according to the directions of the wave vectors (see Table S 13).

S6. Phonon band structures and first Brillouin zones
The following figures show the low-frequency phonon band structures of all studied molecular crystals along with their first Brillouin zones. There, also the high-symmetry points are displayed explicitly to facilitate the recognition of the directions of the corresponding wave vectors with respect to the (reciprocal) lattice.       In addition to all these geometric and material-specific quantities, a dimensionless variable, m, enters Equation (S6), which is only determined by the boundary conditions of the bar, i.e., the kind and position of the supports. 23 As most of the geometric and material-specific parameters in Equation (S6) would not allow a reasonable determination for the given (non-continuous) case of bending molecules, the model is used to estimate frequency ratios rather than absolute frequencies.
Here, we assume that A, , and EI stay constant between the system and that the only parameter that varies is the "beam length", L. The remaining problem is to define the beam length for the bending molecules. In the simplest case of a bending bar, the transverse bending eigenmodes show nodes (i.e., positions of no transverse displacement) at the positions of the supports. That means, that a reasonable way to define the beam length in the molecules is to find those nodes and define the lengths between them as the "effective length", Leff. An analysis of the bending vibrations (at ) reveals that this effective length approximately obeys the following empirically found equation: Here, dring is the (short) diameter of the hexagonal ring of benzene and n is an integer ranging from 1 (benzene, 1A) to 5 (pentacene, 5A) and, thus, counts the number of rings in the molecule.
Additionally, a parameter, s, enters the above equation to account for the observation that the ratio of Leff and the total molecular length decreases for longer acenes. For the bending modes, s=0.5 is used, whereas for the intramolecular in-plane bending modes (for which the same model based on a bending beam is employed), s=0.7 was found suitable to approximate Leff in the atomic displacement patters and, thus, the associated frequencies.
Once Leff is known, the following equation can be used to estimate the frequency ratios relative to a reference system: Here, we used naphthalene (2A) as a reference for the mode, as this molecule is the simplest case which has a defined long molecular axis and shows no phonon hybridisation effects.

S7.2 Intramolecular torsion model
Similarly to the case of bending molecules introduced above, we employed a simple model of a torsional oscillator -with a torsional stiffness, ktor, and the torsional moment of inertia, Itor -to describe the frequency evolution as a function of the number of rings, n: As a next step, the two ingredients, ktot and Itor, should be discussed. To this end, it is vital to understand how the intramolecular torsional mode distort the molecular geometries. In those modes, one can clearly observe a nodal line, which corresponds to the y-axis as shown in Figure S 16. All atoms above the y-axis (i.e. all atoms with z > 0) undergo torsional motion around the zaxis in one direction as one largely rigid unit, while all the atoms below the y-axis (z < 0) undergo a torsion in the other direction (as one largely rigid unit). Depending on whether the number of rings, n, of the molecule is even (n mod 2 = 0) or odd (n mod 2 = 1), different atoms must be excluded from the total moment of inertia around the z-axis, zz (this quantity is introduced in Equation (S13) in Section S7.3): for even n, two central carbon atoms do not contribute to Itor because they lie on the nodal plane and do not move, while for odd n, two carbons and to hydrogens must be excluded (see Figure S 16). For the sake of completeness, the analytic expression for Itor used for the acenes with n rings is given in Equation (S10). (S10) In addition to the total (rotational) moment of inertia around the z-axis, zz (as introduced in Equation (S13) below), Itor depends on the atomic masses of the C and H atoms, mC and mH, and the (average) nearest-neighbour C-C and C-H distance, rCC and rCH. Note that the numeric values (1.396933 Å and 1.089885 Å) for these distances were obtained from the average nearestneighbour distances in benzene (1A) and were kept constant for all estimations of Itor for the sake of simplicity.
Besides this odd-even effect in the calculation of Itor, also the torsional stiffness, ktor, is supposed to depend on the precise location of the nodal plane. Since the atoms with y > 0 and those with y < 0 undergo torsions as largely rigid units, the strongest local bond distortions are observed close to the nodal plane (y = 0). Therefore, we consider the local bond geometries in the proximity of the nodal plane as the restoring "springs" of this vibration. In case of odd-ringed acenes, this directly affects four C-C bonds, which are distorted from the energetically preferred planar (aromatic) arrangement. We consider this number of distorted bonds in the torsional stiffness as four times an arbitrary torsional reference stiffness,  0 tor. In case of even-ringed acenes, also the C-C bond in the centre of the molecule finds itself in an environment of non-planar bonds. Since this bond is "doubly unsatisfied", as the bonding environment on both sides is changing in opposite directions, this bond is counted twice such that ktor equals six times  0 tor in this case. To sum up, the torsional stiffness, ktor, is modelled to depend on the number of rings, n, in the following way: , = 2 0 ⋅ (3 − mod 2) (S11)

S7.3 Rotational rigid intramolecular modes model
In analogy to the analysis of the translational rigid intermolecular modes (TRIMMs) presented in the main text, a similar analysis was applied also to the rotational rigid intermolecular modes (RRIMMs). For a rotational vibration, a suitable measure for the inertia of the oscillator is the appropriate moment of inertia, ii, which must be calculated for rotations around all axes separately (see below). If this quantity is known, the effective rotational molecular stiffness, Krot, can subsequently be calculated from the (angular) frequencies:  (S15) As already discussed in the main text, this calculation of the moments of inertia is somewhat problematic as the axes around which the rotations are observed do not always perfectly agree with the symmetry axes of the molecules, for which the expressions in Equation (S13) to (S15) hold.
Nevertheless, using the numeric values obtained from the above equations still allows to draw several conclusions when analysing the resulting effective rotational force constants as a function of the associated moments of inertia. This analysis is shown in Figure S   The evolution of the rotational force constants for RRIMMs corresponding (roughly) to rotations around the short molecular axes (see Figure S 17(b)) show much more unclear trends. Although for the antiphase RRIMM, a linear function through the origin seems to be suitable to describe the evolution (with antracene and naphthalene lying distinctly below the trend), for the in-phase RRIMM, a first-order polynomial (i.e., slope and offset) were fitted. In this case, a clear trend of decreasing frequency with the molecule length is observed (by a factor of ~1.6 in total from 1A to 5A). Therefore, the rotational force constants are not supposed to depend linearly on the moments of inertia (which would result in a constant frequency). Instead, due to the positive offset, b0, and the slope, b1, the frequencies are expected to decrease with √ 1 + 0 / . However, direct connections to the molecular length are hampered by the fact that in contrast to zz, yy depends on the number of rings per molecule in a non-linear way.
Finally, for the last type of RRIMMs, again a more or less linear relationship between xx and the calculated rotational force constant can be observed (suggesting roughly constant frequencies); the deviations from the linear trends are, however, larger than in all other cases discussed so far.
Still, it seems that the increase in the magnitude of the effective rotational force constants can keep up with the (superlinear) increase in the moment of inertia (approximately) associated with this kind of motion.

S8. Sound velocities (long-wavelength limit)
For the spatial distributions of sound velocities, the group velocities of the acoustic bands were calculated for a small sphere (with radius equalling 1 % of the length of the shortest reciprocal lattice vector in each system to guarantee comparability; see Table S 14) on a 360×180 mesh of azimuthal, , and polar angles, , respectively. For the sake of completeness, the average, maximum, minimum and standard deviation of the TA(1,2) and LA sound velocities on these reciprocal-space spheres are listed in Table S 14, and the deviations from the respective averages are plotted as projections onto unit spheres for all studied systems in the following.
Note that, since the sound velocities , were sampled at meshes in spherical coordinates, average values, 〈 〉, and standard deviations, std( ),were calculated in the following way (with p = TA1/TA2/LA)

S9. Effects of Polymorphism
Here, the most pronounced differences in the phonon-derived quantities introduced above and in the main text will be briefly discussed comparing monoclinic (1A) and orthorhombic benzene (1A-o) as well as the pentacene polymorphs I (5A) and II (5A-II).

S9.1 Orthorhombic benzene
As the orthorhombic phase of benzene (1A-o) contains twice the number of molecules per unit cell, it is impossible to assign the inter-and intramolecular modes on a one-to-one base as it was done for the remaining systems. However, one still can find interesting differences in many phonon-derived quantities discussed above and in the main text.
Regarding the phonon band structure (see Figure S  The a priori probability for avoided crossings is significantly reduced in the orthorhombic polymorph as the system exerts a space group with much more IRREPs. 27 Thus, it is less likely that two bands with the same symmetry meet, which is the prerequisite for phonon hybridisation. In order to emphasise this aspect, Figure S 21 shows the phonon band structure of 1A and 1A-o along one exemplary path with the bands coloured according to their assigned IRREPs. It can be seen that the LA band in 1A meets two bands with the same IRREP along its path and, therefore, experiences two avoided crossings. Conversely, in 1A-o, the LA band only meets one band with the same IRREP along its way to Z so that only one avoided crossing is observed. Also the optical modes ≤ 4 THz (along this path) undergo drastically fewer avoided crossings although the density of bands is notably increased in this region in 1A-o. to the former. This perfect face-to-edge herringbone packing of the molecules (meeting at nearly 90°) in this system seems to be especially suitable to allow a fast propagation of compression waves in this direction.

S9.2 Pentacene polymorph II
The direct comparison is much easier for the two considered polymorphs of pentacene. As the polymorph II (5A-II) shows certain differences in the packing arrangement (see Section S2), especially the intermolecular low-frequency phonons are expected to show some differences to the polymorph I (5A). As far as low-frequency intramolecular bands are concerned, Figure S 22 shows that their frequencies stay essentially the same. The largest for these modes is on the order of 0.3 THz (~4.05 THz in 5A and 3.75 THz in 5A-II) and is observed for the in-phase second-order bending modes, while all the differences for the remaining modes are smaller by at least a factor of ~3. This is most probably a result of the fact that in this mode in 5A, one of the two molecules shows a strong rotational character in addition to the second-order bending due to hybridisation with a closely lying long axis RRIMM. As expected, the changes of the intermolecular modes' frequencies are notably larger, where they amount to up to 0.83 THz for the difference in antiphase RRIMM around the -plane normal.

Figure S 22: Evolution of selected intramolecular and intermolecular frequencies at
Also the in-phase RRIMM around the short axis experiences a significant drop in frequency, which is so large in magnitude (~-0.67 THz) that it comes to lie below the a3/long TRIMM. We attribute this to the much larger geometric offset between two layers of molecule, short, in 5A-II (see Table   S 3). As a result, the molecules have more space to rotate around the short axis and, thus, the interaction strength (effective force constant) must be smaller, resulting in a lower frequency.
Interestingly, the TRIMMs are much less by the different polymorphic packing than the RRIMMs: the two higher ones (a1/normal and a2/short) essentially keep their frequency, while the frequency of the a3/long TRIMM increases by ~0. 16 THz. Analysing the animations, this seems to be the consequence of the increased bending character that the a3/long TRIMM in 5A additionally exhibits compared to the equivalent TRIMM in 5A-II.
Regarding the dispersion of the bands in both polymorphs, the low-frequency band structures (see Figure S 14 and Figure S 15) show many similarities. Even the band widths are typically nearly identical. One of the most pronounced differences apart from the slightly changed order of low-frequency intermolecular bands is the increased band width of the TA modes in 5A-II in X and Y compared to 5A. The larger band widths can also be seen in the low-frequency zoom of the DOS in Figure S 6: here the typical quadratic (Debye-like) onset of the DOS (as a result of the linear dispersion of the acoustic bands) increases slightly more gradually in 5A-II. Apart from that, hardly any (significant) differences in the low-frequency DOS can be perceived, let alone in the higher intramolecular region > 7 THz. As a result, thermodynamic properties at elevated temperatures, for which also intramolecular modes contribute to a significant extent, are expected to exhibit nearly no differences between the two polymorphs (such as the heat capacities discussed in Section S10).
Also regarding the sound velocities, 5A and 5A-II are very similar. The average LA and TA2 sound velocities are essentially identical with even the standard deviations being nearly the same.
A more pronounced difference is observed in the TA1 sound velocities: both the average TA1 sound velocity and the associated standard deviation are slightly higher in 5A than in 5A-II, although the distributions of group velocity differences look practically identical for 5A and 5A-II (see Figure  Concerning the group velocities of the optical modes, the maximum of the low-frequency DOGV is slightly shifted to lower values, implying that among the bands with frequencies ≤ 7 THz, the group velocities are somewhat reduced in 5A-II compared to 5A.

S10. Molar heat capacities
It turns out that at 300 K, all the studied systems display nearly the same value of the normalised heat capacity, CV/(3NkB), (~0.3324) as discussed in the main text. Based on this observation, one can crudely estimate the molar heat capacity, which is typically the one measured in experiments, in the following way: values are always slightly underestimated. One reason for this is that, in general, the relation Cp > CV holds at finite temperatures, since CV implicitly neglects the influence of thermal expansion. 33 In addition to the aspects already discussed in the main text it can be seen that polymorphism does not seem to have a notable effect on the heat capacity at constant volume. The intramolecular modes, which increasingly gain importance at higher temperatures, are too similar so that they obscure the pronounced but comparably few differences in the low-frequency region. Only at very low temperatures (see Figure S 23(c)), one can see that CV/(3NkB) is slightly smaller for 1A-o than for 1A. This observation emphasises the findings discussed in Section S9.1: the weight of the lowfrequency (intermolecular) DOS is slightly shifted to higher frequencies. Therefore, somewhat higher temperatures are needed to thermally activate those phonons' contribution to the heat capacity.
Conversely, the (normalised) heat capacity of the 5A and 5A-II does not even show significant differences in this low-temperature region. Obviously, the few shifts in frequencies and the few differences in band dispersion are too minor to achieve notable variations in CV/(3NkB).  28,29 , naphthalene 29,30 , anthracene 29,31 , tetracene 32 , and pentacene 32 in the corresponding colours. (c) Zoom into the region indicated by the dashed box. Note that the symbols on top of the lines do not correspond to the actually calculated data points (those lie much more densely), but rather serve as guides to the eye.

S11. Comparison of low-frequency intramolecular between the crystals and the isolated molecules
Interesting insight into the non-covalent interactions due to the molecular packing in the studied organic crystals can be obtained by comparing the frequencies of eigenmodes in the isolated molecules to the equivalent -phonons in the crystalline systems. The molecular frequencies in the isolated molecules were calculated using the Gaussian 16 package 34