Coarse-grained lattice dynamics calculations combined with independent stiffness approximation: a comparative study on polymorphic molecular crystals

Abstract

We applied a previously reported lattice dynamics calculation scheme based on molecular vibration analysis, combining a coarse-graining (CG) scheme with an independent stiffness approximation (ISA) method, to polymorphs of p-dichlorobenzene crystals as a test case. To extend the applicability of the theoretical framework to complex crystal structures, we exploited the molecular symmetry within the lattice so that the ISA could effectively reduce the computational cost to a manageable level. We carefully examined the influence of the coupling between phonon modes on the accuracy of their frequencies under a truncated dimension of the CG space. Having obtained the phonon band structures of each of the three polymorphs, their thermodynamic functions were calculated, which demonstrated a subtle but meaningful difference originating from the crystal structure.

---

Introduction

The lattice dynamics of molecular crystals have gained increasing importance in recent crystal engineering studies.¹ With the widespread adoption of density functional theory (DFT), the theoretical calculation of collective lattice vibrations (phonons) has become easier.^2–5 Currently, mainstream computational methods include the finite-displacement supercell approach (frozen-phonon method)^6–8 and density functional perturbation theory (DFPT).^9–12 One of the fundamental differences between them lies in the distinct forms of the force constants obtained during the initial calculations: the force constant matrix, K, in real space versus the dynamic matrices, D, at specific wave vectors in reciprocal space. The sizes of these matrices correspond to the degrees of freedom of the atomic vibrations, endowing these approaches with both high precision and broad applicability. However, for molecular crystals containing many atoms within a unit cell, atomistic-level calculations significantly increase the computational costs. Moreover, these approaches output results as individual atomic motions, requiring a posteriori analysis to relate atomic displacements to the vibrational behavior of the molecule as a consolidated rigid body.¹³ Therefore, by applying this method to molecular crystals, several attempts have been made to reduce the computational cost.^14,15

To address these issues, we developed a method for reducing the dimensions of the atomic displacement vectors of normal-mode vibrations, referred to as the coarse-graining (CG) scheme.^16–21 In this approach, a stiffness matrix, Φ, which describes the degrees of freedom associated with molecular rigid-body motions and vibrations, is readily obtained by transforming the eigenvalue matrix with normal-mode coordinates of a molecular assembly, without explicitly deriving the force constant matrix, K. Similar to the finite displacement method, we need molecular clusters that cover the desired expanse, with which, however, we can accelerate computation using the analytical differentiation protocol as utilized in the DFPT method. Moreover, we recently proposed a method called the independent stiffness approximation (ISA),^22,23 which divides the whole into pieces, thereby converting a multibody computing problem into a multigroup two-body computing problem. Using the CG-ISA scheme, we previously completed phonon calculations for a simple molecular crystal (the β phase of p-dichlorobenzene).^24–26 Through systematic theoretical validation and comparison with other computational results, the CG-ISA scheme was successful for simple systems whose constituent molecules could be related by translational symmetry. However, the CG-ISA approach is additionally challenged by molecular vibrations in more complex symmetric relationships and increased degrees of freedom when each unit cell contains two or more molecules (Z). These challenges are directly reflected in the structure of Φ, which is closely linked to the symmetry of molecular crystals. The following section presents a new theoretical extension of the CG-ISA scheme for systems with Z > 1, as demonstrated by illustrative computational examples. In this paper, we chose p-dichlorobenzene again for a cases study for the following reasons: (1) a simple small molecule yet has three distinct polymorphs, of which crystal structure data at the same temperature are available; (2) basically van der Waals crystals without a complicated networks of hydrogen bond; (3) experimental phonon band structures (in part) and Raman data are available. These features merit our understanding of the structure–property relationship in molecular crystalline materials as an initial stage of our study.

Theory

Outline of coarse-graining method

The vibration of the particles in a lattice can be understood within the harmonic approximation by solving the eigenvalue equation of the mass-weighted Hessian matrix, where the key parameters correspond to the inertia matrix, M, and stiffness matrix, K. In the atom-based analysis, the elements of these matrices correspond to the degrees of freedom of the atomic vibrations.^27–29 For a molecular assembly containing N atoms, the matrix dimensions are 3N × 3N. By introducing a coarse-graining matrix, B, which projects the atomic displacement vectors onto a reduced space spanned by a basis set representing rigid-body (three translational and three rotational) motions (T_x, T_y, T_z, R_x, R_y, R_z) and α intramolecular vibrational modes (V₁,…,V_α), the stiffness matrix is transformed to correspond to the vibrational modes of the molecule with a size of 6 + α per molecule (eqn (1)–(3) and Fig. 1).^16,17 We refer to the number 6 + α as the dimension of the CG basis.

Γ⁻¹ = B^TMB (1)

Φ = B^TKB (2)

B = (T_x, T_y, T_z, R_x, R_y, R_z, V₁, …, V_α) (3)

Fig. 1 Transformation of the stiffness matrix in the coarse-graining scheme. The subscripts of K denotes atoms, whereas subscript of Φ denotes lattice nodes.

The eigenvalue equation equivalent to the equation of motion can be expressed using the matrices U and Ω², which collect the eigenvectors u and ω², respectively, as in eqn (4):

Γ^1/2ΦΓ^1/2U = UΩ² (4)

The elements of the eigenvectors directly correspond to the vibrational modes of the molecule, distinguishing them from traditional methods of using an atom as the unit of vibration. This approach significantly facilitated the analysis of molecular vibrational modes. By leaving out intramolecular vibrational modes from the CG basis, low-frequency vibrational modes can be approximated, which aids in reducing the computational costs for large molecular systems.^19,20 When α is set to 0, the system adopts the minimal basis, approximating the molecule as a rigid body and considering only intermolecular translational and rotational modes. When α is set to 3N − 6, the system adopts the full basis, fully describing the vibrational modes corresponding to 3N degrees of freedom.

For an arbitrary matrix, A, we define notation [A] for its Fourier transformation using eqn (5), where {A}_σ,σ′ denotes a partial matrix corresponding to an off-diagonal block between the constitutional units labeled σ and σ′.²¹

In eqn (5), a_j denotes the j-th lattice vector that serves as the basis for the linear expansion of n₁a₁ + n₂a₂ + n₃a₃ for an arbitrary lattice node, R_σ, of the molecule σ. Furthermore, P denotes the number of molecules included in a single period, and k is assumed to be a continuous reciprocal vector at the limit of P → ∞. Using this expression, we formulated the dynamic equation (eqn (6)) under periodic boundary conditions (PBC).

[Γ^1/2][Φ][Γ^1/2][U(k)] = [U(k)][Ω(k)²] (6)

Independent stiffness approximation

In practical phonon calculations based on the CG scheme, information on all the blocks of the stiffness matrix, Φ, associated with the central lattice node (numbered 0) is required. The PBC is introduced via the Fourier transformation of blocks to obtain [Φ] (eqn (7)).

To completely count the molecules in a period, we defined σ = 9n₁ + 3n₂ + n₃, a decimal expression of a ternary number “n₁n₂n₃”. In this numbering system, the central 1 lattice node and surrounding 26 lattice nodes (n₁, n₂, n₃ = {−1, 0, 1}) are labeled with integers from −13 to 13.²³

To obtain each block {Φ}_0,σ with a reduced computational load, we introduced the concept of ISA, wherein the total intermolecular force acting on a molecule can be calculated as the sum of pairwise interactions, assuming that the stiffness between any two molecules is independent of all others.³⁰ When the number (Z) of molecules contained in each node is 1, the off-diagonal block {Φ}_0,σ can be calculated by selecting two molecules that are related by a translational operation. When Z > 1, a single node contained multiple molecules related to inversion, mirror reflection, rotation, glide, and helical operations. Under the CG space of 6 + α dimensions, each node of the crystal includes Z × (6 + α) vibrational modes, implying that the size of Φ increases accordingly. Considering Z = 2 as an example, the composition of {Φ}_σ,τ is shown in Fig. 2. Here, different molecules within a single node are designated as A and B, which are convertible through a certain symmetry operation allowed within a given space group. By combining the letters A, B, C, … for different operations with node index σ, different molecules within each node can be uniquely identified. For simplicity, we introduce a notation such as {φ}_0A,σB, denoting the (A, B) subblock of the {Φ}_0,σ block (eqn (8)).

{φ}_0X,σY ≡ {{Φ}_0,σ}_X,Y (X,Y = {A,B,C, …}) (8)

Fig. 2 Structures of the stiffness matrix under different definitions of vibrational units: transformation from a lattice node containing two molecules (left) to that containing one molecule (right).

For Z > 1, two approaches to fill out {Φ}_0,σ are possible. One approach is to expand the range of the vibrating unit in the ISA cluster by treating all Z molecules as a single unit (Fig. 2, left side). This treatment results in a “mobile pair” comprising 2Z molecules, simplifying the subsequent theoretical problem to be equivalent to a Z = 1 system. However, as Z increases, the computational cost for the normal-mode calculation increases roughly in proportion to (2Z)³ ∼ (2Z)⁴, and the dimension (α) of the intramolecular vibration mode needs at least 6Z − 6 to include the mutual displacement of Z molecules. Another approach, which was adopted in the present study, is to construct {Φ}_0,σ by “tessellating” smaller {φ}_0X,σY subblocks corresponding to every two molecules within a node (Fig. 2, right-hand side). As Z increases, the number of sub-blocks increases in proportion to Z², and their arrangement becomes increasingly intricate. Because a single ISA cluster calculation yields only one {φ}_0X,σY sub-block, the total number of calculation routines increases proportionally with Z². Minimizing the number of calculation routines while correctly constructing {Φ}_0,σ is the core theoretical challenge in applying the CG-ISA scheme to systems with Z > 1.

In the CG-ISA scheme, where each discrete molecule is assigned to a vibrating unit, σX, the numerical relationship in the {Φ}_0,σ blocks fully reflects the symmetry of a given molecular crystal by setting an appropriate coordinate system. As described in the previous study on a system with Z = 1, the transpose of {Φ}_0,σ is identical to {Φ}_0,−σ: this relation reduces the number of off-diagonal blocks to half. For approximately spherical molecules, one molecule typically has 12–18 adjacent molecules that contribute significantly to its binding energy, meaning that only 6–9 sets of molecular pairs must be calculated. A similar relationship (eqn (9)) also holds between sub-blocks that correspond to arbitrary molecular pairs.

{φ}^T_0X,σY = {φ}_σY,0X = {φ}_0Y,−σX = {φ}^T_−σX,0Y X,Y = {A, B, C, …} (9)

For molecules in the translational relation (X = Y), the relation in eqn (9) can halve the number of sub-blocks necessary for the calculation. For molecules in non-translational relations (X ≠ Y), an appropriate coordination transformation of the {φ}_0X,σY subblock can afford another subblock that corresponds to a molecular pair sharing the equivalent geometrical relationship. This transformation relation was formulated and verified theoretically, the details of which will be presented in a future study. This relation reduces the number of independent subgroups to 1/h, where h is the order of the space group.

To fully exploit eqn (9) in reducing the computational cost, one needs to specify molecular pairs (0X′, σ′Y′) that are in a geometrical relation equivalent to a given molecular pair (0X, σY). For example, for the calculation of the γ phase of p-dichlorobenzene, six ISA clusters shown in Fig. 3 are essential, where molecules 0A and 0B are convertible by the 2₁-helix operation (along the b-axis). For the molecular pair (0A, 3A), namely, σ = 3 (a = 0, b = 1, c = 0) and X = Y = A (identical operation), the C₂ operation only reverses the directions of the a- and c-axes, resulting in a σ′ value of 3 (because both a and c are zero). Thus, we found that the molecular pair (0B, 3B) has a geometrical relation equivalent to (0A, 3A), whereas {φ}_0B,3B is obtained by the coordination transformation of {φ}_0A,3A using P[C₂], a matrix of two-fold rotation around the b-axis (eqn (10)).

{φ}_0B,3B = P[C₂]{φ}_0A,3AP[C₂]^T (10)

Fig. 3 ISA clusters for molecular pairs, as designated by 0X, σY combinations.

With a series of such transformation relations, the force field exerted on only one molecule is sufficient to line up the sub-blocks necessary for all Z molecules in a lattice node. Therefore, an increase in Z does not increase the number of ISA clusters necessary for the calculation.

As listed in Table 1, by computing six ISA clusters, the matrix elements for 30 {φ}_0X,σY subblocks are finally obtained. In this table, the mobile pairs (0X, σY) correspond to the ISA clusters shown in Fig. 3. By calculating the coarse-grained stiffness matrix, subblocks {φ}_0X,σY for the “equivalent” molecular pair are obtained. The transpose of these subblocks corresponds to the molecular pairs in the “transpose” column. The transformation by P[C₂] gives the subblocks for the pairs in the “C₂ rotation” column, and their transpose corresponds to the last column. As the number of symmetry operations in the system increased, the coordination transformation steps became more intricate. Further details are beyond the scope of this study. Hence, transformations within the CG-ISA framework will be discussed in future publications.

Table 1 List of equivalent molecular pairs (0X, σY) of the γ form of p-dichlorobenzene, for which subblocks {φ}_{0X, σY} are calculated from the ISA cluster containing the designated mobile pair. See text for molecular pairs in the other columns

Mobile pair (symmetry operation)	Equivalent	Transpose	C2 rotation	C2 rotation + transpose
0A, 0A (identical)	0A, 0A	—	0B, 0B	—
0A, 3A (b translation)	0A, 3A	0A, −3A	0B, 3B	0B, −3B
0A, 7A (a − b + c translation)	0A, 7A	0A, −7A	0B, −13B	0B, 13B
0A, 10A (a + c translation)	0A, 10A	0A, −10A	0B, −10B	0B, 10B
0A, 0B (21 helix)	0A, 0B	0B, 0A	0B, 3A	0A, −3B
	0B, 1A	0A, −1B	0A, −4B	0B, 4A
0A, 9B (21 helix + a translation)	0A, 9B	0B, −9A	0B, −6A	0A, 6B
	0B, 10A	0A, −10B	0A, −13B	0B, 13A

---

As shown in eqn (8) or (9), our method is formulated without any restriction on the Z number, which means that the scheme applies to various crystals with an arbitrary space group. However, in the actual computation of such a crystal, we must present a full lineup of P matrices necessary for the space group under consideration. Meanwhile, it should be noted that the situation is not easy for a crystal system that has two or more independent units. In that case, we must treat the units separately, which means that the number of possible ISA clusters rapidly increases and the selection of the clusters becomes more complicated. At present, we cannot find an efficient way to select the mobile pairs essential to the phonon band calculation. That is more of a mathematical problem than a chemistry problem.

Calculation

The scheme of the whole procedure was illustrated in Fig. S1 (SI). We selected the α and γ phases of p-dichlorobenzene as calculation targets (Fig. 4). The α phase, CCDC refcode DCLBEN07, monoclinic, P2₁/a, Z = 2, a = 14.664 Å, b = 5.74 Å, c = 3.925 Å, and β = 111.77°. The γ phase, CCDC refcode DCLBEN03, monoclinic, P2₁/c, Z = 2, a = 8.624 Å, b = 6.021 Å, c = 7.414 Å, β = 127.51°.²⁴ For both cases, the inversion point was located at the center of the molecule; hence, only half of the molecules were independent. For computational convenience, we redefined the crystal system as the P2₁ space group such that a lattice node comprises two molecules. With this treatment, only “translation” and “2₁-helix” are left as the symmetry operations of the crystal system. Otherwise, the “inversion” operation would make a molecule overlay with itself, losing the uniqueness of the molecular orientation, hence the uniqueness of the CG matrix. Despite this symmetry reduction, the original symmetry could still be utilized to reduce the number of independent molecular pairs. For example, owing to the centrosymmetry of the molecule, the geometrical relationship of (0A, 0B) is identical to (0B, 1A), which should be distinguished in the P2₁ space group (Fig. S2). We confirmed a substantial equivalence between {φ}_0A,0B and {φ}_0B,1A within a trivial error, which may be attributed to the nonstrict symmetry in the truncated ISA clusters. Table 1 lists the molecular pairs in the “equivalent” column.

Fig. 4 Unit cell structures of p-dichlorobenzene polymorphs: (a) α form (refcode: DCLBEN07) and (b) γ form (DCLBEN03).

Upon constructing the ISA clusters, the essential mobile pairs should be selected based on the distance between the gravity centers of molecules (d_Mole), the interatomic proximity (d_Atom), and the trace of the {φ} subblock, as measures of intermolecular interactions.²³ For the β-DCB, molecules with d_Atom > 5 Å were ignored, and this selection had almost no effect on the phonon band profiles. In the present two systems, where similar intermolecular interactions occur, we used the same criteria to determine the molecules that comprise mobile pairs in the ISA clusters: 0A, 3A, 7A, 10A, 0B, and 9B for γ-DCB (Table 1) and 0A, 1A, 3A, 4A, 0B, and 1B for α-DCB.

The structures of the constructed ISA clusters were iteratively optimized until they reached a substantial minimum potential energy within the given molecular environment.^31,32 The calculations were performed using the Gaussian16 program³³ at the B3LYP-D3/6-31G(d) level of theory, a hybrid DFT augmented with an empirical dispersion correction.³⁴ The same level of theory was applied to the normal-mode calculation for the ISA clusters, in which the mobile pair was set to move freely, and its surrounding molecules were frozen. The coarse-grained stiffness constants were obtained via an appropriate process using normal coordinate data, as described in our previous studies.^16,17 An increase in the Z value of the system does not necessarily imply that the size of the ISA clusters must be expanded. In this study, we employed clusters of sizes comparable to those used in our previous study on DCLBEN06 (Z = 1). Errors originating from the truncation of the cluster size were compensated for by the acoustic sum rule (ASR) correction,^27,28 redefining the translation-relevant part of the diagonal block to maintain consistency with the isotropy and homogeneity of the space.

Fig. 5 shows the first Brillouin zone (FBZ) of DCLBEN03 (γ phase) and DCLBEN07 (α phase), along with the k-paths that connect the high-symmetry points. As understood from the characteristics of reciprocal space, the shapes of the FBZs can differ depending on the lattice constants, even though the two crystalline systems share the same space group. The selection of the k-path in the FBZ shape determines the horizontal axis. Hence, it has a primary influence on the appearance of the phonon band but does not alter the density of states (DOS) and subsequent calculations.

Fig. 5 FBZ and k-paths connecting high-symmetry reciprocal points: (a) α phase and (b) γ phase of p-dichlorobenzene.

Results and discussion

Effect of the CG-space dimension

The {φ}_0X,σY sub-blocks of the stiffness matrix were calculated by normal-mode analysis performed at the B3LYP-D3/6-31G(d) level of theory, which was demonstrated in our previous study²³ to provide sufficiently reliable results for the β phase crystal. Then, {Φ}_0,σ blocks were constructed according to the crystal symmetry and subsequently transformed to a reciprocal space representation [Φ] using Fourier convolution under the assumption of PBC (see Theory section). Fig. S3 shows the low-frequency region of the phonon band (black lines) obtained within the framework of the CG-ISA scheme (minimal basis; α = 0). Fig. 6 compares this result with non-CG (full basis; α = 3N − 6 = 30) calculation (blue lines). A full-range representation of the full-basis results is provided in the SI (Fig. S4). The similarity between the two sets of results demonstrates the feasibility of applying the CG-ISA method to the calculations for systems with Z = 2 and Z = 1. Namely, the tessellation of {φ}_0X,σY and the ASR correction were effective when the single lattice node comprised two molecules. Note that owing to the increased number of molecules within the lattice, the number of phonon branches doubled to 12 up to 140 cm⁻¹, as compared with β-DCB.

Fig. 6 Comparison of phonon bands and DOS between the minimal basis of CG space (α = 0, black) and full basis (α = 3N − 6 = 30, blue). (a) γ phase and (b) α phase.

Close inspection of the γ-DCB result (Fig. 6(a)) revealed a notable discrepancy between minimal- and full-basis results around the #11 and #12 branches (∼120–140 cm⁻¹), while the frequency distribution in α-DCB (Fig. 6(b)) was nearly identical. To assess the quality of the CG space for a mobile pair in a single ISA cluster, we employed the fidelity index (FI),^19,20 which is a measure of the Frobenius inner product of the original and coarse-grained mass-weighted displacement matrices. It is empirically known that the CG space is of satisfactory quality when the FI is 0.8 or higher. For α-DCB and γ-DCB, the average FI values at the minimal-basis condition were 0.93 and 0.98, respectively, both achieving a sufficient level. The observed discrepancy suggests the need for deeper consideration in addition to the reproducibility of Φ.

When combined with ISA, the CG calculation necessarily involves multiple clusters with various FI values, indicating that the approximation quality cannot be determined based solely on this parameter. In other words, the CG space (selection of a limited number of intramolecular vibration modes) suitable for one ISA cluster may not always be appropriate for all other clusters. The previous research on the β-DCB crystal demonstrated the excellent agreement between the minimal-basis and full-basis results for the low-frequency phonon band and DOS. This crystal system may have been an ideal case in which the constituent molecule was sufficiently stiff to separate the frequency regions between the rigid-body motion and intramolecular vibration modes.

To determine the factors affecting the approximation quality of the CG-ISA scheme, we examined the influence of low-frequency intramolecular modes. Fig. 7 shows the phonon bands for modes #1–#16, calculated using the full-basis condition. For γ-DCB, modes #11–14 share the same frequency region, while for α-DCB, those modes are well resolved into {#11, #12} and {#13, #14} sets. In quantum mechanics theory, microstates with the same symmetry interact and change each other's frequency due to the avoided crossing (AC) rule,^35–37 which is accompanied by the mixing of different states. If a selected set of bases in the CG space lacked some modes relevant to the AC effect, mode mixing would not impact the phonon bands, giving rise to a considerable displacement from the full-basis results. This assumption was validated by the calculation of phonon band and DOS with α = 1 (Fig. 8), where we can see an excellent agreement with those under the full-basis condition. This result demonstrates the necessity of including intramolecular vibrational modes that share the same frequency region while maintaining a sufficiently high FI value.

Fig. 7 Phonon bands for modes #1 to #16 in full basis: (a) γ phase and (b) α phase. Modes #11–14 are shown in purple, yellow, red, and gray, respectively.

Fig. 8 Comparison of phonon bands and DOS between CG basis α = 1 (green) and full basis (blue) of γ phase.

Comparison of the Γ-point frequency

The full-basis CG result is the “fine-grained” limit and hence does not reduce computational accuracy. Therefore, a comparison between the reduced- and full-dimensional CG spaces allows for the evaluation of the approximation quality of the CG scheme. Meanwhile, the arbitrariness in constructing ISA clusters influences computational accuracy. In the ISA method, errors arise mainly from the inequality of the molecular environment among clusters owing to the truncation of the range of molecular placement and the subtle disturbance of symmetry owing to structural optimization. In our current computational protocol, the introduction of ASR corrections can compensate for a major proportion of the errors stemming from the aforementioned reasons. However, for systems with more complex space groups that require a larger number of different types of ISA clusters, other unexpected errors may occur. These issues will be discussed in detail elsewhere. At the present stage, we only compared the results of the moderately sized cluster calculations with those of other methods and the experimental values.

Given the lack of reliable experimental or theoretical phonon band diagrams in the literature, we validated our computational results by comparing the frequencies at the Γ-point (Table 2). For comparison with the experimental results, we corrected the frequencies using a scaling factor³⁸ corresponding to the level of calculation. Overall, we observed good agreement between the calculated and experimental values, whereas significant discrepancies were observed in certain modes. Specifically, for α-DCB, the CG approximation performed well, yielding only minor differences between the experimental and calculated results (α = 1 or 30). Among them, the frequencies of modes #6–9 are closer to those of modes #7–10 obtained from the reference calculation. This “misalignment” may have arisen from differences in the definition of mode ordering, necessitating a labelling of the modes using irreducible representations based on group theory. In γ-DCB, the CG-ISA scheme demonstrated better agreement, with experimental values below 100 cm⁻¹, particularly for mode #8, which was much improved as compared to the reference calculations. However, for modes #11 and #12, the calculated values did not correspond to the experimental data. Instead, the calculated frequencies 134 and 141 cm⁻¹ for modes #13 and #14, respectively, are in good agreement with the experimental values for modes #11 and #12.

Table 2 Comparison of frequencies (cm⁻¹) at Γ-point and experimental data, along with the values calculated using another method

Mode#	α form				γ form
	CG-ISA^a		Calcd.^b	Obsd.^c	CG-ISA^a		Calcd.^b	Obsd.^c
	α = 1	α = 30			α = 1	α = 30
a Modes were numbered according to the increasing order of frequency at the Γ-point.b Calculated data from ref. 39 and 40.c Raman data from ref. 41 (T = 1.2K).
1–3	0	0	—	—	0	0	—	—
4	27.3	25.0	35.8	27	48.5	48.8	44.6	—
5	37.0	37.2	29.5	33	51.0	50.3	44.8	—
6	52.3	53.1	45.8	46	52.3	52.2	52.5	52
7	60.9	60.3	53.2	—	64.1	63.1	65	67
8	66.2	65.8	57	59	73.8	73.5	55.2	76
9	69.4	68.1	67	66	83.3	83.1	81.2	86
10	70.9	71.3	66	67	90.8	91.3	102.4	—
11	102.6	101.1	108.5	109	112.9	112.7	114.1	133
12	111.3	110.9	117.5	117	121.0	120.7	111.1	143

---

Changes in the lattice structure primarily affect low-frequency phonons, that is, intermolecular vibrational modes. In contrast to the moderate variation shown in the low-frequency band, the high-frequency phonon branches were nearly identical across the three polymorphs, indicating that the intramolecular vibrational modes in the p-dichlorobenzene crystals were minimally affected by changes in the crystal structure (Fig. 9). This insensitivity is attributed to the absence of strong intermolecular interactions, such as hydrogen bonds, within these crystals, along with the small size and high rigidity of the molecule. Nevertheless, a closer inspection of the differences revealed some common effects of the molecular surroundings in the crystalline environment. Fig. 9 shows a plot of the eigenwave numbers of the polymorphs subtracted by that of the isolated p-dichlorobenzene molecule. The difference ranged from −20 to +40 cm⁻¹, and the wavenumber-dependent trend of the displacement was roughly similar among the three polymorphs. It is notable that in the low-frequency region (up to 300 cm⁻¹), the plots were similar between the α and γ forms, while, in the middle-frequency region (800–1000 cm⁻¹), the plots were similar between the α and β forms.

Fig. 9 Comparison of frequency at Γ-point in the high-frequency region (characterized as intramolecular modes). Given that the β phase has an FBZ half the size of those of the other phases, each frequency appears twice.

Comparative studies on DOS and thermochemistry

Currently, comparative studies of polymorphs are of central importance in crystal engineering.^42–45 In the context of crystal structure prediction, intensive efforts have been made to calculate the relative stability of crystalline polymorphs with high precision, considering both entropic and enthalpic contributions.^46–50 Nevertheless, listing polymorphs in the correct stability order remains difficult. In the following, we assess the applicability of our method to explain the relative stability of polymorphs of molecular crystals. The phonon DOS functions of the three phases (α, β, and γ-DCBs) in the low-frequency region were compared (Fig. 10). The frequency distributions were quite similar among the three phases below 35 cm⁻¹, the range primarily dominated by the acoustic modes. Given the absence of strong intermolecular interactions such as hydrogen bonding, the vibrational modes in this range are hardly affected by changes in the molecular crystal structure.

Fig. 10 Phonon DOS in the low-frequency region (characterized as intermolecular modes) of three polymorphs of p-dichlorobenzene.

In the range of 35 to 100 cm⁻¹, the vibrational modes are primarily low-frequency optical modes. This frequency region for the β phase with Z = 1 is dominated by molecular libration, while the other two phases exhibit a mixing of intramolecular vibration modes, such as out-of-plane C–Cl bending motions. For the α and γ-DCBs, distinct DOS peaks appear near the boundary of the acoustic and optical branches. At these frequencies, phonon modes with considerably low group velocities tend to accumulate. By contrast, the β phase does not exhibit a pronounced DOS peak, indicating a more loosely packed molecular arrangement, which leads to a broader distribution of phonon frequencies. The vibrational modes transition from lattice motion to intramolecular vibrations above 100 cm⁻¹, where the intramolecular branches of the α and γ-DCBs show different DOS peaks, depending on their respective crystal structures.

Experimentally, the order of thermodynamic stability is β > α > γ: the γ form transitions to the α form at 271 K with a transition enthalpy (Δ_trH) of 1.256 kJ mol⁻¹, then transitions to the β form at 304 K with a Δ_trH of 0.214 kJ mol⁻¹ before finally melting at 326 K.^51,52 The transition entropy (Δ_trS) was accordingly calculated to be 4.6 and 0.70 J K⁻¹ mol⁻¹ for the γ-to-α and α-to-β transitions, respectively. One of the primary milestones of this study is explaining this order of stability in terms of thermochemistry. As described above, provided that sufficient dimensions were used, the CG scheme had a negligible influence on the frequency distribution, indicating that its effect can also be negligible on other phonon-derived properties, such as thermodynamic properties. Fig. 11 presents the thermochemical functions of the three phases of p-dichlorobenzene. As established by the statistical dynamics of molecular systems, vibrations contribute to constant-volume heat capacity, and hence entropy and Helmholtz/Gibbs free energy. Although these functions differed only slightly among the three phases, we found that non-negligible differences arose from the low-frequency phonon contributions to the total functions (Fig. 11(a)).

Fig. 11 Thermochemical functions of the three polymorphs calculated with full-basis ISA condition. (a) Heat capacity at constant volume, and (b) entropy, where the insets show the low-frequency contributions with dashed lines. (c) Gibbs free energy relative to the γ form, where the enthalpies of the α and β forms were assumed to be +0.376 and +0.820 kJ mol⁻¹, respectively, relative to the γ form.

The entropy values at ∼100 K were in the order β > α > γ (Fig. 11(b)), indicating that the −TS term in the free energy contributes to stabilizing the β phase as the temperature increases. However, the difference was so small that we could not draw reasonable G–T curves using the experimentally obtained enthalpies (Δ_γ→αH = 1.256 kJ mol⁻¹, Δ_α→βH = 0.214 kJ mol⁻¹): an error of 0.1 kJ mol⁻¹ in ΔH causes a displacement of 40–60 K in T at the intercept of the curves. In a future study, we plan to calculate the changes in enthalpy among the different phases, for example, taking into account the electronic energies of the constituent molecules and intermolecular interaction energies among the molecules, while it was revealed that the calculation requires an accuracy of 0.01 kJ mol⁻¹ to predict transition temperatures reasonably. Instead, we attempted to obtain the ideal enthalpy that yielded the experimentally determined transition temperature (Fig. 11(c)). As a result, Δ_γ→αH and Δ_α→βH were estimated to be 0.376 and 0.820 kJ mol⁻¹, respectively, which were in moderate agreement with the observed values. The sum of the two is 1.196 kJ mol⁻¹, moderately close to the corresponding experimentally determined value of 1.470 kJ mol⁻¹.

For more rigorous quantification, we should consider the temperature-, pressure-, or stress-dependent anisotropic strain of a crystal lattice, which induces changes in the phonon band diagram that are not straightforward to predict.^53–56 Strictly speaking, the thermochemical functions calculated from a single-phonon band diagram do not exactly reproduce their temperature dependence; however, this issue is not specific to the CG-ISA method.^57,58 Given these limitations, the current thermochemistry results restrict the analysis to a qualitative one. Moreover, while the interactions between phonons play a crucial role in various properties, such as thermal conductivity, the harmonic approximation affords “static” phonons, which lack inter-phonon interactions and hence cannot accurately reproduce “dynamic” thermal phenomena. A conventional solution requires incorporating the anharmonicity effects into DFT calculations.^59–64 The present CG-ISA scheme was established based on harmonic approximation, and the modification of the theory to incorporate the anharmonicity effects awaits numerous challenges.

Conclusions

In this study, we enhanced the CG-ISA scheme for lattice dynamics calculations in molecular crystal systems with Z > 1. By analyzing the relationship between the subblock components {φ}_0X,σY of the stiffness matrix and crystal symmetry, a comprehensive theory was proposed to cover the construction of the ISA clusters and tessellation of Φ. The current progress has enabled the calculation of the α and γ phase crystals of p-dichlorobenzene with Z = 2. We evaluated the approximation quality of the CG scheme by comparing the calculations for different dimensions of the CG space. By selecting a sufficient number of intramolecular vibrational modes considering the avoided crossing rule, the truncated basis result was successfully reproduced. The calculated frequencies at the Γ-point are in good agreement with the experimental values, while the assignment of phonon branches based on group theory is necessary for rigorous matching with them. Comparisons of the three polymorphs revealed that the acoustic and high-frequency intramolecular modes were minimally affected by the lattice structure, whereas the low-frequency optical modes exhibited notable sensitivity to the lattice structure. Similar vibrational modes in different crystalline phases exhibit frequency shifts owing to structural variations. Thermochemical calculations demonstrated that the three polymorphs resulted in only a subtle difference but were qualitatively meaningful, encouraging us to perform more accurate evaluations by considering additional effects.

Author contributions

All the authors contributed to the conception and design of the study. Yue Wang mainly conducted the computation, data processing, and data analyses. Hirohiko Houjou organized the whole research project. The two collaborated on the theoretical content of the study. Y.W. prepared the first draft, and H.H. reviewed and edited the manuscript. All the authors approved the final version of the manuscript.

Conflicts of interest

The authors declare no competing financial interest.

Data availability

Data for this article, including Gaussian input files, output files, and coarse-graining analysis files are available at OSF at https://osf.io/583ux/overview?view_only=fc04744fdbd74e62af4f60c2a8fab9dd.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d6cp00159a.

Acknowledgements

This work was partially supported by JST SPRING, Grant Number JPMJSP2108, and World-leading Innovative Graduate Study program for Materials Research, Information, and Technology (MERIT-WINGS) at the University of Tokyo.

Notes and references

1. L. Catalano, K. M. Hutchins, C. J. Bardeen and M. T. Ruggiero, Lattice Dynamics: The Unexplored Multidimensional Dynamic Playground of Molecular Crystalline Materials, Cryst. Growth Des., 2024, 24, 2301–2303.
2. G. M. Day, J. A. Zeitler, W. Jones, T. Rades and P. F. Taday, Understanding the Influence of Polymorphism on Phonon Spectra: Lattice Dynamics Calculations and Terahertz Spectroscopy of Carbamazepine, J. Phys. Chem. B, 2006, 110, 447–456.
3. M. R. C. Williams, D. J. Aschaffenburg, B. K. Ofori-Okai and C. Schmuttenmaer, Intermolecular Vibrations in Hydrophobic Amino Acid Crystals: Experiments and Calculations, J. Phys. Chem. B, 2013, 117, 10444–10461.
4. F. Zhang, H.-W. Wang, K. Tominaga and M. Hayashi, Mixing of intermolecular and intramolecular vibrations in optical phonon modes: terahertz spectroscopy and solid-state density functional theory, WIREs Comput. Mol. Sci., 2016, 6, 386–409.
5. F. Zhang, H.-W. Wang, K. Tominaga and M. Hayashi, Intramolecular Vibrations in Low-Frequency Normal Modes of Amino Acids: L-Alanine in the Neat Solid State, J. Phys. Chem. A, 2015, 119, 3008–3022.
6. A. Togo and I. Tanaka, First Principles Phonon Calculations in Materials Science, Scr. Mater., 2015, 108, 1–5.
7. A. Togo, L. Chaput, I. Tanaka and G. Hug, Phys. Rev. B:Condens. Matter Mater. Phys., 2010, 81, 174301.
8. K. Parlinski, Z. Q. Li and Y. Kawazoe, First-Principles Determination of the Soft Mode in Cubic ZrO2, Phys. Rev. Lett., 1997, 78, 4063–4066.
9. S. Baroni, S. de Gironcoli, A. Dal Corso and P. Giannozzi, Phonons and Related Crystal Properties from Density-Functional Perturbation Theory, Rev. Mod. Phys., 2001, 73, 515–562.
10. X. Gonze, First-Principles Responses of Solids to Atomic Displacements and Homogeneous Electric Fields: Implementation of a Conjugate-Gradient Algorithm, Phys. Rev. B:Condens. Matter Mater. Phys., 1997, 55, 10337–10354.
11. S. Baroni, P. Giannozzi and A. Testa, Green's-Function Approach to Linear Response in Solids, Phys. Rev. Lett., 1987, 58, 1861–1864.
12. P. Giannozzi, S. de Gironcoli, P. Pavone and S. Baroni, Ab Initio Calculation of Phonon Dispersions in Semiconductors, Phys. Rev. B:Condens. Matter Mater. Phys., 1991, 43, 7231–7242.
13. L. Legenstein, L. Reicht, T. Kamencek and E. Zojer, Anisotropic Phonon Bands in H-Bonded Molecular Crystals: The Instructive Case of α-Quinacridone, ACS Mater. Au, 2023, 3, 371–385.
14. T. Kamencek, S. Wieser, H. Kojima, N. Bedoya-Martínez, J. P. Dürholt, R. Schmid and E. Zojer, Evaluating Computational Shortcuts in Supercell-Based Phonon Calculations of Molecular Crystals: The Instructive Case of Naphthalene, J. Chem. Theory Comput., 2020, 16, 2716–2735.
15. L. Soprani, A. Giunchi, M. Bardini, Q. N. Meier and G. D’Avino, Accurate and Efficient Phonon Calculations in Molecular Crystals via Minimal Molecular Displacements, J. Chem. Theory Comput., 2025, 21, 8073–8085.
16. H. Houjou, Coarse Graining of Intermolecular Vibrations by a Karhunen-Loève Transformation of Atomic Displacement Vectors, J. Chem. Theory Comput., 2009, 5, 1814–1821.
17. H. Houjou, Evaluation of coupling terms between intra- and intermolecular vibrations in coarse-grained normal-mode analysis: Does a stronger acid make a stiffer hydrogen bond?, J. Chem. Phys., 2011, 135, 154111.
18. H. Houjou, Modelling intra- and intermolecular vibrations under the harmonic oscillator approximation: from symmetry-adapted to coarse-grained coordinate approaches, J. Math. Chem., 2017, 55, 532–551.
19. M. Isogai and H. Houjou, Indices to evaluate the reliability of coarse-grained representations of mixed inter/intramolecular vibrations, J. Mol. Model., 2018, 24, 221.
20. M. Isogai, M. Seshimo and H. Houjou, Optimizing a coarse-grained space for approximate normal-mode vibrations of molecular heterodimers, J. Mol. Model., 2021, 27, 140.
21. H. Houjou and M. Seshimo, Formulation of a Phonon Band Calculation for Molecular Crystals Using a Coarse-Grained Coordinate Approach under Periodic Boundary Conditions, J. Math. Chem., 2022, 60, 613–636.
22. H. Houjou, Y. Wang and S. Okamura, Formulation of Phonon Band Calculation Scheme Using Intermolecular Stiffness Matrix Represented upon Coarse-grained Coordinate System, J. Comput. Chem. Jpn., 2021, 20, 147–149.
23. Y. Wang, M. Seshimo and H. Houjou, Size-Reduction of Phonon Band Calculation for Coarse-Grained Molecular Crystals Using “Independent Stiffness Approximation”, J. Chem. Theory Comput., 2025, 21, 6048–6060.
24. G. L. Wheeler and S. D. Colson, Intermolecular interactions in polymorphic p-dichlorobenzene crystals: The α, β, and γ phases at 100 K, J. Chem. Phys., 1976, 65, 1227–1235.
25. M.-M. Thiéry and C. Rérat, Calculation of crystal and molecular structures of the temperature and pressure polymorphs of para-dichlorobenzene p-C₆H₄Cl₂, J. Chem. Phys., 2003, 118, 11100–11110.
26. P. A. Reynolds, J. K. Kjems and J. W. White, Lattice vibrations in chlorobenzenes: Experimental dispersion curves for β-paradichlorobenzene by neutron scattering, J. Chem. Phys., 1974, 60, 824–834.
27. S. Califano, V. Schettino and N. Neto, Lattice Dynamics of Molecular Crystals (Lecture Notes in Chemistry, 26), Springer, 1981.
28. M. T. Dove, Introduction to Lattice Dynamics, University Press, 1993.
29. E. B. Wilson Jr., J. C. Decius and P. C. Cross, Molecular Vibrations, Dover Publications, 1980.
30. V. Coropceanu, R. S. Sánchez-Carrera, P. Paramonov, G. M. Day and J.-L. Brédas, Interaction of Charge Carriers with Lattice Vibrations in Organic Molecular Semiconductors: Naphthalene as a Case Study, J. Phys. Chem. C, 2009, 113, 4679–4686.
31. R. Koibuchi, I. Yoshikawa, Y. Shigemitsu and H. Houjou, Comprehensive Computational Analysis of the Polymorph-Dependent Chromism of N-Salicylideneaniline Using a Cluster Model of Molecular Crystals, Cryst. Growth Des., 2024, 24, 238–251.
32. R. Koibuchi, M. Makita, H. Shingai, I. Yoshikawa and H. Houjou, Controlling the stability of tautomeric polymorphs in a multistep proton-transfer system using a homologue approach to fine-tune the potential energy landscape, Cryst. Growth Des., 2025, 25, 7655–7668.
33. M. J. Frisch, G. W. Trucks, H. B. Schlegel, et al., Gaussian 09w (Revison D), Gaussian Inc, Wallingford, 2019.
34. S. Grimme, S. Ehrlich and L. Goerigk, Effect of the Damping Function in Dispersion Corrected Density Functional Theory, J. Comput. Chem., 2011, 32, 1456–1465.
35. M. Baggioli, B. Cui and A. Zaccone, Theory of the phonon spectrum in host-guest crystalline solids with avoided crossing, Phys. Rev. B, 2019, 100, 220201.
36. J. Baumert, C. Gutt, V. P. Shpakov, J. S. Tse, M. Krisch, M. Müller, H. Requardt, D. D. Klug, S. Janssen and W. Press, Lattice dynamics of methane and xenon hydrate: Observation of symmetry-avoided crossing by experiment and theory, Phys. Rev. B:Condens. Matter Mater. Phys., 2003, 68, 174301.
37. Y. Zhu, B. Wei, J. Liu, N. Z. Koocher, Y. Li, L. Hu, W. He, G. Deng, W. Xu, X. Wang, J. M. Rondinelli, L.-D. Zhao, G. J. Snyder and J. Hong, Physical insights on the low lattice thermal conductivity of AgInSe₂, Mater. Today Phys., 2021, 19, 100428.
38. CCCBDB listing of precalculated vibrational scaling factors. https://cccbdb.nist.gov/vibscalex.asp (accessed 12–25, 2025).
39. H. Bonadeo, E. D’Alessio, E. Halac and E. Burgos, Lattice dynamics, thermodynamic functions, and phase transitions of p-dichloro- and 1,2,4,5-tetrachlorobenzene, J. Chem. Phys., 1978, 68, 4714–4721.
40. J. Serrier, F. Brehat, B. Wynckeet and A. Hadni, Vibrations de réseau de quelques dérivés dihalogénés du benzène, Can. J. Chem., 1979, 57, 1814–1822.
41. A. P. J. M. Jongelenis, T. H. M. van der Berg, J. Schmidt and A. van der Avoird, Vibron band structure in chlorinated benzene crystals: lattice dynamics calculations and Raman spectra of 1,4-dichlorobenzene, J. Phys.: Condens. Matter, 1989, 1, 5051.
42. L. Yu, Polymorphism in Molecular Solids: An Extraordinary System of Red, Orange, and Yellow Crystals, Acc. Chem. Res., 2010, 43, 1257–1266.
43. J. P. Reddy, D. Swain and V. R. Pedireddi, Polymorphism and Phase Transformation Behavior of Solid Forms of 4-Amino-3,5-dinitrobenzamide, Cryst. Growth Des., 2014, 14, 5064–5071.
44. A. Zimmermann, B. Frøstrup and A. D. Bond, Polymorphs of Pridopidine Hydrochloride, Cryst. Growth Des., 2012, 12, 2961–2968.
45. Q. Zhang and X. Mei, Two New Polymorphs of Huperzine A Obtained from Different Dehydration Processes of One Monohydrate, Cryst. Growth Des., 2016, 16, 3535–3542.
46. C. Greenwell, J. L. McKinley, P. Zhang, Q. Zeng, G. Sun, B. Li, S. Wen and G. J. O. Beran, Overcoming the difficulties of predicting conformational polymorph energetics in molecular crystals via correlated wavefunction methods, Chem. Sci., 2020, 11, 2200–2214.
47. A. J. Zaczek, L. Catalano, P. Naumov and T. M. Korter, Mapping the polymorphic transformation gateway vibration in crystalline 1,2,4,5-tetrabromobenzene, Chem. Sci., 2019, 10, 1332–1341.
48. A. M. Reilly and A. Tkatchenko, Role of Dispersion Interactions in the Polymorphism and Entropic Stabilization of the Aspirin Crystal, Phys. Rev. Lett., 2014, 113, 055701.
49. J. Hoja, A. M. Reilly and A. Tkatchenko, First-principles modeling of molecular crystals: structures and stabilities, temperature and pressure, WIREs Comput. Mol. Sci., 2017, 7, e1294.
50. J. Hoja, H.-Y. Ko, M. A. Neumann, R. Car, R. A. DiStasio Jr. and A. Tkatchenko, Reliable and practical computational description of molecular crystal polymorphs, Sci. Adv., 2019, 11, eaau3338.
51. A. Dworkin, P. Figuiere, M. Ghelfenstein and H. Szwarc, Heat capacities, enthalpies of transition, and thermodynamic properties of the three solid phases of p-dichlorobenzene from 20 to 330 K, J. Chem. Thermodynam., 1976, 8, 835–844.
52. H. A. J. Oonk, A. C. G. van Genderen, J. G. Blok and P. R. van der Linde, Vapour pressures of crystalline and liquid 1,4-dibromo- and 1,4-dichlorobenzene; lattice energies of 1,4-dihalobenzenes, Phys. Chem. Chem. Phys., 2000, 2, 5614–5618.
53. O. Hellman, I. A. Abrikosov and S. I. Simak, Temperature dependent effective potential method for accurate free energy calculations of solids, Phys. Rev. B:Condens. Matter Mater. Phys., 2013, 87, 104111.
54. I. Errea, M. Calandra and F. Mauri, Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: Application to platinum and palladium hydrides, Phys. Rev. B:Condens. Matter Mater. Phys., 2014, 89, 064302.
55. A. Wang, J. Yin, F. A. Goudreault, M. Côté, O. Hellman and S. Poncé, Antagonistic impact of thermal expansion and phonon anharmonicity on the phonon-limited resistivity of elemental metals from first principles, arXiv, preprint, arXiv:2511.17478 DOI:10.48550/arXiv.2511.17478.
56. M. T. Ruggiero, S. Ciuchi, S. Fratini and G. D’Avino, Electronic Structure, Electron-Phonon Coupling, and Charge Transport in Crystalline Rubrene Under Mechanical Strain, J. Phys. Chem. C, 2019, 123, 15897–15907.
57. E. J. Blancas, Á. Lobato, F. Izquierdo-Ruiz, A. M. Márquez, J. M. Recio, P. Nath, J. J. Plata and A. Otero-de-la-Roza, Thermodynamics of solids including anharmonicity through quasiparticle theory, npj Comput. Mater., 2024, 10, 267.
58. D. B. Zhang, T. Sun and R. M. Wentzcovitch, Phonon Quasiparticles and Anharmonic Free Energy in Complex Systems, Phys. Rev. Lett., 2014, 112, 058501.
59. M. Hutereau, P. A. Banks, B. Slater, J. A. Zeitler, A. D. Bond and M. T. Ruggiero, Resolving Anharmonic Lattice Dynamics in Molecular Crystals with X-Ray Diffraction and Terahertz Spectroscopy, Phys. Rev. Lett., 2020, 125, 103001.
60. B. Monserrat, N. D. Drummond and R. J. Needs, Anharmonic vibrational properties in periodic systems: energy, electron-phonon coupling, and stress, Phys. Rev. B:Condens. Matter Mater. Phys., 2013, 87, 144302.
61. J. C. A. Prentice and R. J. Needs, Using forces to accelerate first-principles anharmonic vibrational calculations, Phys. Rev. Mater., 2017, 1, 023801.
62. D. Mitoli, A. Erba, V. Barone and M. Mendolicchio, Anharmonicity in Molecular Crystals: Generalized Perturbation Theory Meets Periodic Computations, J. Phys. Chem. Lett., 2025, 16, 9956–9962.
63. A. Erba, J. Maul, M. Ferrabone, P. Carbonnière, M. Rérat and R. Dovesi, Anharmonic Vibrational States of Solids from DFT Calculations. Part I: Description of the Potential Energy Surface, J. Chem. Theory Comput., 2019, 15, 3755–3765.
64. A. Erba, J. Maul, M. Ferrabone, R. Dovesi, M. Rérat and P. Carbonnière, Anharmonic Vibrational States of Solids from DFT Calculations. Part II: Implementation of the VSCF and VCI Methods, J. Chem. Theory Comput., 2019, 15, 3766–3777.

---