Converging ab initio phonon simulations for organic molecular crystals: the effect of charge density grids and phonon dispersion sampling

Abstract

We here explore how some frequently overlooked computational parameters affect the simulation of phonon frequencies in organic molecular crystals within the framework of density functional perturbation theory in a pseudo-core plane wave basis set. Specifically, we investigate how the density of the Fourier grid that is used to map real-space charge density affects the phonon frequencies and eigenvectors. We find that varying the density of this Fourier grid can affect low-frequency phonons by tens of wavenumbers and significantly alter the associated normal mode eigenvectors. Furthermore, we demonstrate that poorly converged charge density representations can lead to substantial errors in simulated thermodynamic quantities, with vibrational free energies affected by 3–4 kJ mol⁻¹ in certain systems. We show how this variation in predicted free energies can have a significant impact on our ability to correctly predict the relative stability of a series of model polymorphic systems. We finally discuss how careful convergence with respect to the Brillouin zone (q-point) sampling is imperative for the correct modelling of phonon dispersion relations in organic molecular crystals, particularly for systems characterised by weak, anisotropic interactions. Whilst no definitive ‘rules of thumb’ emerge for the convergence of these parameters, our findings highlight the critical role they play in obtaining reliable phonon frequencies from density functional perturbation theory. Our results also offer insight into the potential magnitude of errors that could arise in phonon simulations of organic molecular crystals if these parameters are not chosen carefully.

---

1 Introduction

Organic molecular crystals (OMCs) have wide-ranging technological applications, from pharmaceuticals¹ and high-energy-density materials,² to components in micro-electronic,³ micro-optical,⁴ piezoelectric energy-harvesting,⁵ and mechanically actuating devices.⁶ The ability to design, and ultimately optimise, the performance of OMCs for these various technologies requires a clear link between the desired functional behaviour and the underlying molecular and crystal structure. Our typical understanding of ‘structure’ is rooted in representations of static nuclear geometries.⁷ These representations are, of course, extremely useful and underpin most of the last century's developments in the chemical sciences.

As research targets increasingly complex physical phenomena and requires a more detailed characterisation of crystal behaviour, traditional models based on static nuclei are no longer sufficient. For example, to correctly describe the crystal structure landscapes associated with crystal form (polymorph) screening, we must consider the vibrational contributions to the free energy,⁸ which are essential to resolve the subtle energy differences between polymorphs. Similarly, the thermal response of OMCs is also governed by the vibrational behaviour of molecules within the crystalline structure.^9,10 In fact, even seemingly ‘simple’ properties of OMCs, which are often studied by using static-nuclei models, can be affected by vibrational dynamics. For example, a marked renormalisation of electronic band gaps has been observed in molecular crystals once nuclear motion is taken into account.¹¹ Beyond the intrinsic physical properties of crystals, the dynamics of OMCs are central to understanding the chemical reactivity of molecules in the solid state. Mechanisms for the mechanochemical initiation of explosives, for example, have been linked to the channelling of kinetic energy through the vibrational modes of OMCs.^12,13 In reality, most functional properties of OMCs depend on crystal dynamics to some (sometimes subtle) degree. Hence, accurate theoretical models of this behaviour are essential for studying the intricate link between crystal dynamics and OMC properties.

Given the increasingly recognised importance of vibrational dynamics in the study of OMCs, it is not surprising that recent years have seen a significant increase in efforts to simulate this behaviour. Within quantum chemical codes, there are two conventional ways of calculating vibrational dynamics: (1) the ‘brute force’ finite difference method,¹⁴ and (2) the perturbative linear response method.¹⁵ In the former, the force constant matrix is constructed explicitly by computing the forces associated with small atomic displacements in a model crystal structure. This method is relatively straightforward to implement, and is compatible with any energy calculator (including via classical potentials), which has made it particularly popular for modelling phonons. However, because the force constant matrix is calculated explicitly, one can only obtain dynamical matrices (i.e. reciprocal space representations of the force constants) at wavevectors that are commensurate with the number of primitive representations in the supercell that is used in the model. Correspondingly, to generate accurate phonon dispersion curves, large supercells are typically required, quickly rendering the method computationally intractable for large, low-symmetry OMCs. Reducing the cost of finite difference methods for OMC phonon simulations using various approximations or lower-cost force calculators has been the subject of active research.^16,17 In contrast, the linear response method obtains the force constants through a first-order (linear) perturbation of the calculated charge density, and is therefore only compatible with quantum chemical methods. However, the phase of the perturbation can be arbitrarily chosen and hence dynamical matrices at any wavevector are directly available from the primitive unit cell. Therefore, this method is often preferred for simulating the phonon spectra of OMCs.

‘Rules of thumb’ have been explored as strategies for achieving convergence of various calculated properties of OMCs, for example with respect to electronic structure sampling.¹⁸ However, apart from efforts to determine how different levels of theory (e.g., density functionals, force fields, etc.) affect internuclear force constants and hence phonon frequencies,^17,19 little work is available for understanding how calculation parameters themselves can influence the convergence of phonon spectra in OMCs using linear response methods. Given the rapidly growing community of both modellers and experimentalists interested in accessing phonon properties through computational methods, we feel that some discussion of the convergence behaviour, especially of less commonly discussed computational parameters, is timely. Such an investigation is therefore the aim of the present manuscript. We note that it is not our intention to provide a ‘one-size-fits-all’ process to converge phonon frequency calculations, but rather to demonstrate the influence that some less obvious simulation parameters can have on phonon frequencies and common related properties. It is our hope that the present work will provide food-for-thought guidance for those interested in simulating phonon spectra for OMCs, and will help consolidate strategies for obtaining reliable and reproducible phonon data for OMC property prediction.

2 Computational methods

2.1 Input crystal structures

All structures used in this manuscript were obtained from the Inorganic Crystal Structure Database (ICSD)—1,1-diamino-2,2-dinitroethylene (α-polymorph code: 172530, β-polymorph code: 172534)²⁰—or the Cambridge Structural Database: trioxane (REF: TROXAN11),²¹ 1,3,5-triamino-2,4,6-trinitrobenzene (REF: TATNBZ03),²² paracetamol form I (REF: HXACAN36),²³ paracetamol form II (REF: HXACAN08),²⁴ γ-glycine (REF: GLYCIN15),²⁵ α-glycine (REF: GLYCIN04),²⁶ guanidinium nitrate (REF: LETGIA),²⁷ benzoquinone (REF: BNZQUI03),²⁸ benzene (REF: BENZEN20),²⁹ and cubane (REF: CUBANE01).³⁰

2.2 Structure relaxation

Experimental structures were relaxed within the framework of plane-wave density functional theory (DFT) as implemented in CASTEP v22.11.³¹ The Hamiltonian was approximated with the generalised gradient approximation (GGA) functional of Perdew–Burke–Ernzerhof (PBE)³² with the addition of the Grimme D2 dispersion correction.³³ This method has been shown to yield good agreement between calculated and experimental inelastic neutron scattering spectra for CHNO molecular crystals,³⁴ and hence is a reliable approach for modelling phonon spectra in this class of crystals. The nuclear Coulomb potential was attenuated with norm-conserving pseudopotentials as obtained on-the-fly within the CASTEP suite. The wavefunction was expanded in a plane-wave basis set to a kinetic energy cut-off of 1200 eV and sampled on a Monkhorst–Pack k-point grid with 0.05 Å⁻¹ spacing; these parameters were found to be more than sufficient to obtain well-converged total electronic energies and forces for all the model systems. The convergence criteria for the self-consistent field cycles were set for the electronic energy change at <1 × 10⁻¹⁰ eV and for the electronic eigenvalue change at <1 × 10⁻¹² eV. The Broyden density-mixing scheme was used to accelerate the convergence of the SCF calculation, with a mixing amplitude of 0.5 and a maximum mixing G-vector of 1.5 Å⁻¹. Structure relaxation was performed using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) algorithm³⁵ and convergence was accepted once the total energy change reached <2 × 10⁻⁶ eV per atom, residual forces <1 × 10⁻⁴ eV Å⁻¹, ionic displacements <1 × 10⁻⁵ Å and cell strain <0.1 GPa.

2.3 Phonon calculations

All phonon calculations were performed within the linear response method, as implemented in CASTEP v22.11.³¹ In the plane wave linear response formalism, we can define the forces on atom I at position R_I as³⁶where, V_eff is the effective potential exerted on the nucleus from the electron density distribution, F^ion–ion_I describes the force associated with the inter-nuclear interactions, and is the electron density constructed from wavefunctions Ψ_k(r) evaluated at selected wavevectors k in the first Brillouin zone. These wavefunctions are expanded as a Fourier series,with c_k+G representing the coefficients associated with plane wave contributions to Ψ_k(r). The maximum frequency plane wave included in the expansion, defined by the cut-off energy for the basis set E_cut, is limited by . Since ρ(r) is obtained as the square of the wavefunction, its exact representation contains components up to twice the G_max used in the expansion of the wavefunction. Hence, the real-space representation of ρ(r) is defined by a Fourier grid with up to double the point density as compared to the wavefunction grid. In CASTEP, the density of the ρ(r) Fourier grid (the ‘standard grid’) is characterised by G^std_max, where the grid scale relative to the wavefunction grid can be modified by the user. In other solid-state DFT codes, G^std_max is either determined automatically from the user-specified E_cut (Vienna ab initio simulation package,³⁷ VASP) or calculated from a separate kinetic energy cut-off (ABINIT,³⁸ QuantumESPRESSO³⁹).

The charge density stored on a grid defined by G^std_max = 2G_max is exact, and hence provides an exact energy. However, the finite basis set size means it inherently omits the higher-frequency wavefunction components that correspond to the atomic core regions when the wavefunction is pseudised, the representation of which would require a finer sampling of ρ(r). In practice, this can lead to numerical errors when evaluating quantities that depend on the charge density and its spatial derivatives, such as the forces in eqn (1) and, correspondingly, the quality of calculated phonon frequencies. To mitigate this issue, many codes employ a finer grid scale, typically constructed by adding additional G-vectors beyond the formal (2G_max) limit required to represent the charge density. In CASTEP, the scale of this ‘fine’ grid relative to G_max defined by the basis set cut-off energy can be directly set by the user, resulting in the maximum Fourier component G^fine_max = fine grid scale × G_max.

In the present work, we performed linear response phonon calculations using a standard grid defined by G^std_max = 2G_max (corresponding to 35.4943 Å⁻¹ at the selected kinetic energy cut-off of 1200 eV), varying the quality of the fine grid, as described in Table 1. Since the numerical evaluation of the atomic forces depends critically on the sampling of ρ(r), this allows us to investigate the impact of the fine grid scale on the phonon frequencies and eigenvectors calculated for OMCs.

Table 1 Quality of the fine Fourier grid used in the phonon calculations, defined by scale of the fine grid G^fine_max relative to the basis set G_max = 17.7472 Å⁻¹ at the selected E_cut = 1200 eV

Label	Fine grid scale	Size of G^finemax/Å^−1
I	2.0	35.4943
II	3.0	53.2415
III	4.0	70.9887
IV	5.0	88.7359
V	6.0	106.4830
VI	7.0	124.2302

---

To calculate the phonon dispersion relations for our OMCs, their respective Brillouin zones were sampled on Γ-centred Monkhorst–Pack grids with increasingly finer q-point spacing, as discussed in the main text. The phonon frequencies were interpolated from these explicitly calculated grids onto a path between high-symmetry Brillouin zone points, as generated by the SeeKpath utility.^40,41 To study how phonon mode frequencies converge across the dispersion relations, the calculated frequencies were re-ordered at each q-point in order to maximise the magnitude of the scalar product between pairs of eigenvectors that were calculated using different sampling grids. This procedure ensured that our comparisons were made against equivalent eigenvectors (unless otherwise stated in the text).

Vibrational free energies were obtained by performing a standard thermodynamic calculation in CASTEP. Local mode force constants were obtained using an in-house modified version of the solid-state implementation of the LModeA software,^42,43 capable of an automated analysis⁴⁴ of dynamical matrices obtained from Γ-point phonon calculations in CASTEP.

2.4 Simulated inelastic neutron scattering spectra

Experimental inelastic neutron scattering (INS) spectra were obtained from the ISIS INS database for α-1,1-diamino-2,2-dinitroethylene,⁴⁵ trioxane,⁴⁶ 1,3,5-triamino-2,4,6-trinitro-benzene,³⁴ γ-glycine,⁴⁷ α-glycine,⁴⁷ and benzene.⁴⁸ We compared this experimental data to simulated INS spectra obtained using the Mantid 6.8.0 software,^49,50 computed without consideration of higher order quantum events (section S2). All simulations were performed using a temperature of 20 K to best reflect the experimental conditions against which we are comparing.

3 Results and discussion

It is our primary aim to gauge how the construction of the Fourier grid that is used for the real-space charge density (G^fine_max) affects the behaviour of the phonon frequencies and eigenvectors obtained through linear response theory. Our study assumes that the other convergence criteria are well converged. In this respect, we have not performed a complete analysis of how density functionals or dispersion corrections affect the accuracy of the frequencies, though some such analyses are available elsewhere,^17,19 and a brief discussion is given with respect to INS spectra in SI section S2. Moreover, we note that the exact convergence behaviour will depend on the specific pseudopotential construction that one uses. However, we anticipate that our analyses will serve as a useful guide when performing and assessing the reliability of phonon calculations of organic molecular crystals (OMCs).

As model systems for our study, we selected a series of CHNO compounds that reflect common types of bonding in OMCs, Fig. 1. This includes a dispersion-bound crystal (cubane), π-stacked systems (benzene and 1,3,5-triamino-2,4,6-trinitrobenzene (TATB)), examples of crystals with intra-(TATB) and inter-(polymorphs of 1,1-diamino-2,2-dinitroethylene (DADNE), glycine, and paracetamol) molecular hydrogen bonding, and charged systems (glycine and guanidinium nitrate). Within the study, we have also ensured the inclusion of both small, rigid molecules and large, flexible molecules, thereby allowing us to assess how both external and soft internal vibrations respond to the increasing sampling of the charge density.

Fig. 1 Crystal packing diagrams for the eight OMCs studied in this work. In each case, atoms are coloured as oxygen (red), nitrogen (blue), carbon (black), and hydrogen (white). The unit cells are shown as black wires in each case. The coloured shapes alongside OMC names are used throughout the text to identify them, and correspond to CHNO (circle), CHO (square), and CH (triangle) systems. Pairs of polymorphic forms are shown in coloured boxes for DADNE, paracetamol and glycine.

3.1 Phonon convergence at the Brillouin zone centre

Prior to performing phonon calculations, we relaxed all crystal geometries at the PBE-D2 level of theory, subject to the convergence criteria noted in Section 2.2, using a plane wave basis set expanded to 1200 eV (G_max = 17.7472 Å⁻¹). For all cases, we used a standard Fourier grid defined by G^std_max = 2G_max = 35.4943 Å⁻¹, and varied the fine grid G^fine_max as described in Table 1. Using these parameters, the unit cells simulated for each system reproduced the low temperature experimental geometries well, SI section S1. We note that increasing G^fine_max had a negligible effect on the crystal geometries; any change in the phonon frequencies is purely the result of changes in numerical accuracy of the real-space charge densities. For the optimised OMC geometries obtained with the highest fine grid scale (G^fine_max = 124.2302 Å⁻¹, VI in Table 1), we compared simulated inelastic neutron scattering (INS) spectra to experimental spectra available in the ISIS INS Database. This comparison of INS spectra showed good agreement, thereby validating the models used in this study, see SI section S2.

We next explore the convergence of Brillouin zone-centre frequencies by calculating the change in frequency, Δω, as G^fine_max is varied according to Table 1. We compare the frequencies (1) directly in ascending order, as conventionally reported by simulation codes, Fig. 2A (top) and SI section S3.1, and (2) re-ordered to maximise the magnitude of the scalar product between pairs of eigenvectors, to ensure that comparison is made between equivalent normal modes, Fig. 2A (bottom), Fig. 2B and SI section S3.1. As a general observation, we find that only the low-frequency phonon modes are significantly affected by the change in G^fine_max, with negligible effect seen above ca. 400 cm⁻¹ in all cases. In fact, above ω ≈ 400 cm⁻¹, the frequencies change by <1 cm⁻¹ as G^fine_max is increased. This is consistent with the relation ω ∝ , where k is an effective force constant. Within linear response theory, k depends on the response of the charge density ρ(r) to an infinitesimal perturbation λ, k ∝ ∂ρ(r)/∂λ. Consequently, lower frequency phonons, characterised by lower magnitude force constants, are more susceptible to numerical noise in ρ(r).

Fig. 2 Convergence behaviour of Γ-point phonon frequencies for OMCs with increasing G^fine_max. (A) Convergence behaviour of α-DADNE (top) based on a frequency-ordered list and (bottom) using a re-ordered list according to eigenvector. (B) Convergence behaviour of (top) intramolecular hydrogen-bonded crystal TATB and (bottom) dispersion-bound crystal cubane, both using eigenvector re-ordered lists. Plots for all model OMCs are given in SI section S3.

As a representative example, we consider the hydrogen-bonded crystal α-DADNE. When ω are compared in ascending order, an increase in G^fine_max from setting I (G^fine_max = G^std_max = 2G_max) to setting II (G^fine_max = 3G_max) yields Δω with magnitudes up to ±10 cm⁻¹ for the low frequency modes, Fig. 2A (top). The magnitude of Δω decreases with a continued increase in G^fine_max, though frequency variations of >5 cm⁻¹ are still observed up to grid setting IV (G^fine_max = 5G_max). The effect that changing G^fine_max has on the phonon frequencies is magnified when ω are re-ordered according to their corresponding eigenvectors, Fig. 2A (bottom), with a change in G^fine_max from settingI to II yielding Δω > 20 cm⁻¹ in some cases. However, this variation drops again to ca. 2–3 cm⁻¹ once a grid setting of IV is reached.

Fig. 2B compares the influence that G^fine_max has on the phonon frequencies in two different molecular crystals: TATB, which is hydrogen-bonded and contains CHNO ions, and cubane, which is composed only of CH and is stabilised by dispersion forces. As in the case of α-DADNE, Fig. 2A, the lowest phonon frequencies in TATB are highly sensitive to changes in G^fine_max. When normal modes are reordered by their eigenvectors, we observe Δω of up to ca. 30 cm⁻¹ as the grid setting increases from I to II. However, convergence within 2–3 cm⁻¹ is achieved once grid setting IV is used. In contrast, the phonon frequencies of cubane remain virtually unchanged across all grid settings, indicating minimal sensitivity to G^fine_max. This stark difference between the sensitivity of hydrogen-bonded, CHNO-containing systems like DADNE and TATB, and the insensitivity of dispersion-bound CH-containing systems like cubane (and benzene, SI section S3) raises an important question: is this sensitivity of ω to G^fine_max driven only by the specific ions that are involved, or is it also influenced by the nature of the intermolecular interactions?

To delve into this question, we have identified which value of G^fine_max is needed to converge the phonon frequencies in our eight test OMCs to Δω < 1 cm⁻¹, Fig. 3 (top). Equivalent convergence behaviour analyses are shown in SI section S3.2, using convergence thresholds of 0.5, 2 and 5 cm⁻¹. At the outset, we anticipated that higher values of G^fine_max would be needed to converge phonon frequencies for systems that contain harder pseudised cores (here, O ions have the hardest cores). Indeed, we do find that the phonon frequencies for our CH-containing systems (cubane and benzene) are fully converged to within our 1 cm⁻¹ threshold, even with grid setting I. In contrast, much higher G^fine_max are needed to converge the phonon frequencies for CHO (benzoquinone, trioxane) and CHNO (glycine, paracetamol, TATB and DADNE) systems. Of these six OMCs, a grid setting of IV is needed to converge CHNO systems γ-glycine and paracetamol, while a setting of V is needed to achieve convergence for the remaining four systems, albeit with α-DADNE still showing one unconverged low-frequency (ca. 27 cm⁻¹) mode. We do not see any obvious correlation between the nature of intermolecular bonding and the phonon convergence behaviour (note CHO systems that are primarily dispersion-bound also require high G^fine_max to achieve convergence), beyond the general observation that the lowest frequency vibrations are consistently the last to achieve convergence. Moreover, as it is equally challenging to achieve convergence for the rigid molecules trioxane and benzoquinone as for the more flexible molecules DADNE and paracetamol, our findings suggest molecular flexibility is also not a good indicator for the ease of phonon frequency convergence.

Fig. 3 Effect of G^fine_max on the convergence behaviour of phonon simulations for model OMCs. (Top) Summary of convergence behaviour for the OMCs studied in this work, where convergence is considered when Δω < 1 cm⁻¹ as G^fine_max is increased according to the values in Table 1. Unconverged frequencies are shown in white, with converged frequencies displayed in colour. (Bottom) Absolute dot product matrix between phonon eigenvectors calculated at the fine grid setting I (vertical) and fine grid setting VI (horizontal), shown for all phonon modes with ω < 200 cm⁻¹. Tick marks represent mode indices, counting from the bottom-left of each matrix plot.

The change in ω with G^fine_max indicates a change in the strength of the interatomic forces derived from variations in the charge density. Unsurprisingly, we also find that notable changes can arise in the phonon eigenvectors associated with modes of ω < 200 cm⁻¹, Fig. 3 (bottom). This sensitivity of the phonon eigenvectors to G^fine_max is particularly pronounced for the more strongly-bonded OMCs, DADNE, TATB, glycine, and paracetamol. In contrast, weakly bonded crystals cubane, benzoquinone, and trioxane contain only a small number of eigenvectors (two, two, and four, respectively) that are strongly affected by G^fine_max. We do not find that any of the phonon eigenvectors in benzene are affected by changing G^fine_max. These results highlight that also the interpretation of vibrational eigenvectors, whether to establish structure–property relations when interpreting vibrational spectra,⁵¹ to map phonon-driven transformations,^52,53 or otherwise, requires the careful convergence of phonon calculations.

3.2 Effect of Γ-point phonon frequency convergence on physico-chemical properties

We have so far shown that, to obtain converged phonon frequencies and eigenvectors for CHNO-containing molecular crystals in the linear response formalism, one requires considerably finer charge density representations than those dictated by the basis set kinetic energy cut-off. This, and the necessity to perform additional convergence testing, inevitably comes at the expense of increased computational cost. One cannot help but ask the question: is it worth the effort? The answer to this question naturally depends on the purpose of the phonon simulations and the material properties of interest. In this regard, we here explore the effects that convergence of the phonon frequencies with respect to G^fine_max has on phonon-dependent properties. First, we probe the effect of G^fine_max on the vibrational contributions to the crystal free energy, which depend only on the phonon density of states, and hence allows us to examine the importance of an overall converged set of frequencies. Second, we perform local mode analysis (a method to project normal modes onto internal coordinates to obtain effective local mode force constants),⁵⁴ which depends on both the phonon frequencies and the corresponding eigenvectors.

3.2.1 Phonon convergence and vibrational free energy. While density functional theory calculations are inherently performed at 0 K, we can capture some of the free energy's temperature-dependence by including the vibrational free energy using calculated phonon frequencies. Specifically, the vibrational free energy, F_vib, is given bywhere g(ω) is the phonon density of states. The accuracy of these vibrational contributions to the crystal free energy is especially important, given that the energy differences between polymorphic forms are typically on the order of only ca. 2–3 kJ mol⁻¹.^55–57

Using the phonon frequencies obtained in section 2.3, we calculated for each of our OMCs how F_vib change as G^fine_max is increased, Fig. 4 and SI section S4 (note only small changes in zero-point energy with changes in G^fine_max). Stemming from the minimal influence that G^fine_max has on the phonon frequencies for CH crystals, we see that F_vib for cubane and benzene is only negligibly affected by G^fine_max. For these CH systems, increasing G^fine_max from setting I to II was accompanied by a change in F_vib of <0.1 kJ mol⁻¹. In contrast, F_vib for the CHO and CHNO systems is much more sensitive to G^fine_max, again reflecting the sensitivity of their phonon frequencies to this fine grid parameter. In fact, for benzoquinone and TATB, F_vib changes by >2 kJ mol⁻¹—i.e., an amount that is typical of the energy difference between polymorphic forms⁵⁵—when G^fine_max is increased from setting I to II, and 3–4 kJ mol⁻¹ when comparing setting I to VI. For all systems, increasing G^fine_max beyond a setting of III has minimal further effect. Our results demonstrate that errors in F_vib due to insufficient convergence of phonon frequencies with respect to G^fine_max lead to notable change in the overall predicted thermodynamic stability of OMCs, Fig. 4 (middle), and that in some cases this error may be sufficient to cause the re-ordering of predicted polymorph stabilities.

Fig. 4 (Top) Convergence behaviour of the vibrational free energy for OMCs calculated using phonon frequencies obtained with increasing G^fine_max (Table 1). (Middle) Change in the vibrational free energy for OMCs as G^fine_max is increased from setting I (35.4943 Å⁻¹) to setting VI (124.2302 Å⁻¹). (Bottom) Difference in the vibrational free energies calculated for pairs of polymorphs of DADNE, paracetamol and glycine calculated using phonon frequencies obtained with increasing G^fine_max (Table 1). All thermodynamic data are calculated at 300 K.

As a demonstration of this, we compute F_vib for a series of polymorphic forms at a range of G^fine_max settings. We note that the internal energy for our OMCs is unaffected by G^fine_max, and hence the change in the predicted relative polymorph stability becomes dependent only on the relative values of F_vib as G^fine_max is varied, which we plot in Fig. 4 (bottom) and SI section S4. For our model system (glycine) whose F_vib is only minimally affected by G^fine_max, the relative stability of the α- and γ-polymorphs is largely independent of G^fine_max. In contrast, for both DADNE and paracetamol polymorphs, a low, unconverged value of G^fine_max estimates the difference in polymorph stability to be 1–2 kJ mol⁻¹ as compared with predictions done at higher values of G^fine_max.

Using the example of the vibrational free energy, we highlight that bulk thermodynamic properties show similar convergence behaviour with respect to the quality of the charge density as the individual phonon frequencies, despite depending on the overall set of frequencies. In line with our overarching goal to highlight the sensitivity of linear-response phonon calculations to computational parameters, we emphasise the importance of rigorous convergence testing to ensure the reliability of property predictions derived from phonon calculations. It is worth noting that site symmetry splitting, prominent in low symmetry crystals, leads to the presence of a larger number of symmetry independent modes as compared with high symmetry crystals. As such, for low symmetry systems, one might expect that error in the predicted frequencies will have a larger impact on the overall thermodynamic integration.

3.2.2 Local mode analysis. Having seen that G^fine_max can significantly impact on phonon properties that depend on the overall density of states, we now consider a property based on the projection of specific frequencies and eigenvectors: bond strengths, as determined by the effective local mode force constants derived from local vibrational mode theory.⁵⁴ Since local vibrational mode analysis is a bonding descriptor based entirely on vibrational frequencies and eigenvectors as inputs, it serves as an excellent benchmark against which to study the effect G^fine_max has on solid state bonding analysis.

Fortunately, as we have observed in section 3.1, phonons with ω > 400 cm⁻¹ are negligibly affected by the choice of G^fine_max. We might therefore expect that the local vibrational mode analysis of chemical bonding in OMCs will be equally insensitive to G^fine_max. Indeed, local mode analysis of our model OMCs does show only minor variation in the calculated local mode force constants (see SI section S5). For systems found in section 3.1 to require higher G^fine_max values to achieve convergence, we see variation in local mode force constants of no more than ca. 0.1–0.5 mdyn Å⁻¹. This variation is on the order of the differences seen for the local mode force constants associated with the aromatic C–C bonds in different polycyclic aromatic hydrocarbons,⁵⁴ and hence can have impact on bonding studies where fine differentiation between bonding environments is required. However, the local mode force constants were unaffected by G^fine_max for the OMCs whose phonon frequencies also exhibited minimal sensitivity to G^fine_max, such as benzene.

While we find the fortunate situation that local vibrational mode descriptors of chemical bonding are largely insensitive to G^fine_max, this insensitivity is primarily the result of current local mode analysis being restricted to intramolecular bonding analysis. We anticipate the same insensitivity will not hold if local vibrational mode analysis becomes extended to allow for the study of solid state intermolecular interactions in the future.

3.3 Phonon dispersion convergence with respect to Brillouin zone sampling

Having explored the convergence behaviour of Brillouin zone-centre phonon frequencies with respect to the quality of the charge density representation, we now turn our attention to the convergence of the phonon dispersion relations. Phonon dispersion relations describe the dependence of phonon frequencies, ω, on wavevector, q, obtained by diagonalising the dynamical matrix at each sampled q-point in the Brillouin zone. While in theory a complete representation of the phonon dispersion would necessitate an infinite number of dynamical matrices to be calculated, in practice we explicitly calculate the dynamical matrices for a small subset of q-points, and Fourier-interpolate between these points to generate a finer grid. While a sufficiently large G^fine_max ensures well-converged frequencies at individual wavevectors, the phonon dispersion depends on the selection of q-points for which dynamical matrices are explicitly calculated. A sufficiently dense q-point mesh is necessary to accurately interpolate phonon branches across the Brillouin zone. Insufficient sampling can lead to artificial distortions in the dispersion curves, affecting the accuracy of derived properties—but the extent of this effect is not well explored. Here, we examine how the choice of q-point sampling for the explicit computation of the dynamical matrices influences the construction of the phonon dispersion, using a subset of our model OMCs: benzoquinone, γ-glycine, benzene, and cubane. The explicit computation of dynamical matrices was done using G^fine_max values for which Γ-point frequency convergence within Δω ca. <1 cm⁻¹ was achieved (benzoquinone and γ-glycine: V, benzene: II, cubane: I).

Dynamical matrices were explicitly calculated on increasingly finer Γ-centred Monkhorst–Pack (MP) q-point grids, generated based on a gradually decreasing q-point spacing, SI section S6.1, and the force constants were interpolated onto high-symmetry Brillouin zone paths with 0.05 Å⁻¹ step size to produce the phonon dispersion plots, SI section S6.2. For all of our selected OMCs (except cubane), imaginary frequencies were obtained when the phonon dispersion was interpolated from Γ-point phonons only. We found that all phonon branches in our OMCs had positive frequencies when at least two q-points were explicitly sampled along the largest reciprocal space lattice vectors (achieved with MP grid spacing of 0.14 Å⁻¹). This q-point grid spacing (0.14 Å⁻¹) was therefore used as the starting point for our further comparisons.

As we increase the density of q-point sampling from a spacing of 0.14 Å⁻¹ to 0.04 Å⁻¹ (the finest MP grid we sampled), we see substantial variation in the interpolated phonon branch frequencies in the weakly, dispersion-bonded OMCs, cubane and benzene, Fig. 5. This is most prominent in benzene, where variations of Δω > 20 cm⁻¹ are observed across a significant portion of the Brillouin zone, and for branches with average frequencies up to ca. 1100 cm⁻¹. The impact of increasing the q-point sampling density in cubane is somewhat less significant, with Δω on the order of ca. 5–10 cm⁻¹ for branches with average frequency ca. <200 cm⁻¹, but seen across the entire Brillouin zone. In stark contrast, we find that changing the q-point spacing for explicitly calculated dynamical matrices has much less impact on the interpolated phonon branches for benzoquinone and γ-glycine, which contain hydrogen bonds (albeit weak CH⋯O bonds in benzoquinone). In fact, while we do see Δω on the order of 5–10 cm⁻¹ as we decrease the q-point sampling spacing, this remains localised to very small segments of the phonon dispersion curve; across the majority of the dispersion curves for benzoquinone and γ-glycine we only observe Δω on the order of 1–2 cm⁻¹, SI section S6.2. In all the OMCs, it is important to note that the response of the phonon branch frequencies is not linear with increasing q-point grid, with branch frequencies showing much more significant responses to the change in q-point sampling densities for coarser grids. We tend to find that using q-point grids with spacings ≈<0.06 Å⁻¹ yield overall little additional effect on the branch frequencies, albeit with the more weakly bound systems (benzene and cubane) still remaining less well-converged compared to benzoquinone and γ-glycine, SI section S6.

Fig. 5 Change in phonon dispersion as q-point sampling in the Brillouin zone is increased. Phonon frequencies averaged between the coarser and finer q-point spacing are shown as black dots, with the red ‘error bars’ representing the maximum and minimum value of a given frequency at the corresponding q-point, with frequencies matched based eigenvector. Black lines are drawn as a guide for the eye and do not strictly conserve band connectivity. Only selected portions of the phonon dispersion are displayed to represent the convergence behaviour in the low, intermediate, and high frequency range. Further detail is given in SI section S6.2.

The difference in convergence behaviour of the phonon dispersion relation seen across the OMCs can be attributed to the fact that in strongly bonded crystals, the dynamical matrix is well-defined by the presence of several strongly interacting moieties, such as those involved in hydrogen bonding. These strong interactions influence molecular motion over long ranges, leading to phonon dispersions that tend to be less sensitive to changes in the q-point sampling. In contrast, for OMCs where weak interactions play a dominant role, phonon spectra rely on accurately resolving all relevant interatomic force constants, and as a result, greater frequency variations are observed when the Brillouin zone grid spacing is reduced. Whilst we have only a limited dataset, we do not see any obvious correlation between the ease of phonon dispersion convergence and the crystal symmetry. This is reflected in the fact that our low symmetry (monoclinic) system, benzoquinone, converges similarly to our higher symmetry systems, benzene (orthorhombic) and cubane (rhombohedral).

Finally, we turn our attention once again to the vibrational contributions to the free energy of the OMCs, and investigate the variation in F_vib calculated from the stable dispersion relations interpolated from the q-point grids discussed above, Fig. 6. In contrast to the large ΔF_vib observed when Γ-point phonon frequencies calculated with different G^fine_max were used, the largest ΔF_vib in the subset, seen for benzene, is only 0.291 kJ mol⁻¹ despite the notable Δω variations in the dispersion plots. This can be attributed to the phonon frequencies being well-defined at the MP grid points due to a sufficiently fine charge density representation, ensuring that any differences in the interpolation between the explicitly sampled q-points approximately average out. While this is expected for bulk properties like F_vib, more significant errors can be expected in the prediction of anisotropic properties that depend intimately on the degree of dispersion in the phonon branches, like thermal transport,⁵⁸ though this is beyond the scope of the present work.

Fig. 6 Change in the vibrational free energy for OMCs as q-point sampling in the Brillouin zone is increased from the coarsest MP grid yielding a stable dispersion (ca. 0.14 Å⁻¹ q-point spacing) and the finest MP grid used in this study (ca. 0.04 Å⁻¹ q-point spacing), limited by computational cost.

4 Conclusions

We have here explored the convergence behaviour of simulated phonon frequencies and eigenvectors, as obtained through linear response theory in the plane wave density functional theory formalism. In particular, we have probed the sensitivity of both frequencies and eigenvectors to the construction of the charge density grid, a parameter that is often overlooked in literature simulation of phonon properties for organic molecular crystals. We demonstrate that errors on the order of tens of wavenumbers can arise for low-frequency phonons when coarse charge-density grids are used. These errors have significant repercussions for the prediction of thermodynamic quantities. For some of our modelled systems, the errors in thermodynamic quantities that arise from the insufficient convergence of calculated phonon frequencies are on the order of 3–4 kJ mol⁻¹, which is comparable to differences in polymorph stabilities in organic molecular crystals. We demonstrate how the predicted relative stability between polymorphic forms can differ significantly, depending on the degree of convergence in the charge density grid. This finding exemplifies the importance of ensuring sufficient convergence of phonon frequency calculations for studying even ‘simple’ properties of organic molecular crystals. It is worth stressing that to compare polymorphic stabilities one must also take care that the electronic structure descriptors (e.g. k-point sampling and plane wave constructions) are chosen to ensure that the quality of the representation is equivalent for differently sized unit cells. We finally highlight the impact of q-point sampling density on the calculation of phonon dispersion relations for organic molecular crystals, especially in systems with weak, anisotropic intermolecular interactions.

Our studies have been restricted to a single density functional (the generalised gradient approximation functional, PBE) with Grimme's D2 dispersion correction. As the selection of functional and dispersion correction has usually only a small influence on the charge density distribution, and often limited to the low oscillatory components of the charge density structure, we do not anticipate that the convergence behaviour observed in this work will be significantly different for different levels of theory. However, we do stress that the convergence behaviour is intimately linked to the choice of pseudopotential, and hence our work shows only the importance of converging this parameter carefully for a given simulation set-up. Moreover, our work has only considered the impact of this convergence behaviour on the harmonic force constants, but will undoubtedly have an equally important impact on anharmonic force constants. Not only will the inclusion of anharmonic effects influence the accuracy of computed frequencies, but the accurate convergence of anharmonic force constants will be essential should more complex behaviour be sought, including accurate simulations of thermal conductivity and spectral line widths.

Our investigation does not hone in on clear rules for achieving convergence with respect to the charge density representation or q-point sampling. However, we hope that our study will serve as an important demonstration of the significance that these parameters have in the simulation of phonon frequencies within the linear response formalism when using a plane wave basis set. Moreover, we hope that our results contribute to a clearer understanding of the potential magnitude of errors that may be inherent to simulated phonon frequencies where full convergence testing of the demonstrated parameters may be prohibitively expensive.

Conflicts of interest

There are no conflicts to declare.

Data availability

All data associated with this study are available in the supplementary information (SI), or can be obtained at the Zenodo repository, at DOI: https://doi.org/10.5281/zenodo.16631359. Supplementary information is available. See DOI: https://doi.org/10.1039/d5ce01090j.

Acknowledgements

The computations described in this paper were performed using the University of Birmingham's BlueBEAR HPC service, which provides a High Performance Computing service to the University's research community. See http://www.birmingham.ac.uk/bear for more details. Additional calculations were made possible via the UK Materials and Molecular Modelling Hub, which is partially funded by the EPSRC (EP/T022213/1). MM thanks the University of Birmingham for the award of a PhD studentship. TMP acknowledges support from the EPSRC Centre for Doctoral Training in Topological Design (EP/S02297X/1). AALM acknowledges support from a BAM Wilhelm-Ostwald Fellowship.

Notes and references

1. R. Hilfiker and M. v. Raumer, Polymorphism in the pharmaceutical industry: solid form and drug development, Wiley-VCH Verlag, 2nd edn, 2019.
2. S. R. Kennedy and C. R. Pulham, Monographs in Supramolecular Chemistry, Royal Society of Chemistry, 2018, pp. 231–266.
3. R. Samanta, S. Das, S. Mondal, T. Alkhidir, S. Mohamed, S. P. Senanayak and C. M. Reddy, Chem. Sci., 2023, 14, 1363–1371.
4. R. Chinnasamy, J. Ravi, V. Vinay Pradeep, D. Manoharan, F. Emmerling, B. Bhattacharya, S. Ghosh and R. Chandrasekar, Chem. – Eur. J., 2022, 28, e202200905.
5. S. Guerin, S. A. M. Tofail and D. Thompson, NPG Asia Mater., 2019, 11, 10.
6. J. Mahmoud Halabi, E. Ahmed, S. Sofela and P. Naumov, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2020604118.
7. R. Hoffmann and P. Laszlo, Angew. Chem., Int. Ed. Engl., 1991, 30, 1–16.
8. J. Hoja, H.-Y. Ko, M. A. Neumann, R. Car, R. A. DiStasio and A. Tkatchenko, Sci. Adv., 2019, 5, eaau3338.
9. N. Juneja, J. L. Hastings, W. B. Stoll, W. W. Brennessel, S. Zarrella, P. Sornberger, L. Catalano, T. M. Korter and M. T. Ruggiero, Chem. Commun., 2024, 60, 12169–12172.
10. N. F. Xavier, A. M. da Silva and G. F. Bauerfeldt, Cryst. Growth Des., 2020, 20, 4695–4706.
11. A. M. Alvertis and E. A. Engel, Phys. Rev. B, 2022, 105, L180301.
12. J. Bernstein, J. Chem. Phys., 2018, 148, 084502.
13. A. A. L. Michalchuk, J. Hemingway and C. A. Morrison, J. Chem. Phys., 2021, 154, 064105.
14. K. Parlinski, Z. Q. Li and Y. Kawazoe, Phys. Rev. Lett., 1997, 78, 4063–4066.
15. S. Baroni, P. Giannozzi and A. Testa, Phys. Rev. Lett., 1987, 58, 1861–1864.
16. T. Kamencek, S. Wieser, H. Kojima, N. Bedoya-Martínez, J. P. Dürholt, R. Schmid and E. Zojer, J. Chem. Theory Comput., 2020, 16, 2716–2735.
17. J. A. Weatherby, A. F. Rumson, A. J. A. Price, A. Otero De La Roza and E. R. Johnson, J. Chem. Phys., 2022, 156, 114108.
18. J. Binns, M. R. Healy, S. Parsons and C. A. Morrison, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2014, 70, 259–267.
19. F. Brown-Altvater, T. Rangel and J. B. Neaton, Phys. Rev. B, 2016, 93, 195206.
20. J. Evers, T. M. Klapötke, P. Mayer, G. Oehlinger and J. Welch, Inorg. Chem., 2006, 45, 4996–5007.
21. V. Busetti, A. Del Pra and M. Mammi, Acta Crystallogr., Sect. B, 1969, 25, 1191–1194.
22. Z. Chua, C. G. Gianopoulos, B. Zarychta, E. A. Zhurova, V. V. Zhurov and A. A. Pinkerton, Cryst. Growth Des., 2017, 17, 5200–5207.
23. D. A. Druzhbin, T. N. Drebushchak, V. S. Minkov and E. V. Boldyreva, J. Struct. Chem., 2015, 56, 317–323.
24. G. Nichols and C. S. Frampton, J. Pharm. Sci., 1998, 87, 684–693.
25. A. Kvick, W. M. Canning, T. F. Koetzle and G. J. B. Williams, Acta Crystallogr., Sect. B, 1980, 36, 115–120.
26. J. Almlöf, Å. Kvick and J. O. Thomas, J. Chem. Phys., 1973, 59, 3901–3906.
27. A. Katrusiak and M. Szafrański, Acta Crystallogr., Sect. C:Cryst. Struct. Commun., 1994, 50, 1161–1163.
28. M. Bolte, G. Margraf and W. Lerner, CSD Communication, 2002.
29. M. Woińska, S. Grabowsky, P. M. Dominiak, K. Woźniak and D. Jayatilaka, Sci. Adv., 2016, 2, e1600192.
30. R. J. Doedens, P. E. Eaton and E. B. Fleischer, Eur. J. Chem., 2017, 2017, 2627–2630.
31. S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. I. J. Probert, K. Refson and M. C. Payne, Z. Kristallogr. - Cryst. Mater., 2005, 220, 567–570.
32. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
33. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
34. A. A. L. Michalchuk, M. Trestman, S. Rudić, P. Portius, P. T. Fincham, C. R. Pulham and C. A. Morrison, J. Mater. Chem. A, 2019, 7, 19539–19553.
35. R. H. Byrd, J. Nocedal and R. B. Schnabel, Math. Program., 1994, 63, 129–156.
36. S. Baroni, S. de Gironcoli, A. Dal Corso and P. Giannozzi, Rev. Mod. Phys., 2001, 73, 515–562.
37. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
38. X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, P. Ghosez, J.-Y. Raty and D. Allan, Comput. Mater. Sci., 2002, 25, 478–492.
39. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
40. Y. Hinuma, G. Pizzi, Y. Kumagai, F. Oba and I. Tanaka, Comput. Mater. Sci., 2017, 128, 140–184.
41. A. Togo, K. Shinohara and I. Tanaka, Sci. Technol. Adv. Mater.: Methods, 2024, 4, 2384822.
42. Y. Tao, W. Zou, S. Nanayakkara and E. Kraka, J. Chem. Theory Comput., 2022, 18, 1821–1837.
43. F. Bodo, A. Erba, E. Kraka and R. T. Moura, J. Comput. Chem., 2024, 45, 1130–1142.
44. R. T. Moura, M. Quintano, J. J. Antonio, M. Freindorf and E. Kraka, J. Phys. Chem. A, 2022, 126, 9313–9331.
45. A. A. L. Michalchuk, S. Rudić, C. R. Pulham and C. A. Morrison, Chem. Commun., 2021, 57, 11213–11216.
46. B. S. Hudson, Inelastic Neutron Scattering spectrum of 1,3,5-Trioxane, (CH2O)3, measured on the TOSCA instrument, 2017,  DOI:10.5286/edata/668, Accessed Nov 5, 2024.
47. S. A. Rivera, D. G. Allis and B. S. Hudson, Cryst. Growth Des., 2008, 8, 3905–3907.
48. H. Jobic, Stud. Surf. Sci. Catal., 1996, 764, 559–566.
49. O. Arnold, J. Bilheux, J. Borreguero, A. Buts, S. Campbell, L. Chapon, M. Doucet, N. Draper, R. F. Leal, M. Gigg, V. Lynch, A. Markvardsen, D. Mikkelson, R. Mikkelson, R. Miller, K. Palmen, P. Parker, G. Passos, T. Perring, P. Peterson, S. Ren, M. Reuter, A. Savici, J. Taylor, R. Taylor, R. Tolchenov, W. Zhou and J. Zikovsky, Nucl. Instrum. Methods Phys. Res., Sect. A, 2014, 764, 156–166.
50. Mantid Project, Mantid 6.8.0: Manipulation and Analysis Toolkit for Instrument Data, 2024,  DOI:10.5286/SOFTWARE/MANTID6.8, Accessed: March 11, 2025.
51. S. Hunter, P. L. Coster, A. J. Davidson, D. I. A. Millar, S. F. Parker, W. G. Marshall, R. I. Smith, C. A. Morrison and C. R. Pulham, J. Phys. Chem. C, 2015, 119, 2322–2334.
52. X. Liu, A. A. L. Michalchuk, B. Bhattacharya, N. Yasuda, F. Emmerling and C. R. Pulham, Nat. Commun., 2021, 12, 3871.
53. A. J. Zaczek, L. Catalano, P. Naumov and T. M. Korter, Chem. Sci., 2019, 10, 1332–1341.
54. E. Kraka, W. Zou and Y. Tao, WIREs Comput. Mol. Sci., 2020, 10, e1480.
55. J. Nyman and G. M. Day, CrystEngComm, 2015, 17, 5154–5165.
56. J. Nyman, L. Yu and S. M. Reutzel-Edens, CrystEngComm, 2019, 21, 2080–2088.
57. M. Arhangelskis, F. Topić, P. Hindle, R. Tran, A. J. Morris, D. Cinčić and T. Friščić, Chem. Commun., 2020, 56, 8293–8296.
58. H. Yang, J.-Y. Yang, C. N. Savory, J. M. Skelton, B. J. Morgan, D. O. Scanlon and A. Walsh, J. Phys. Chem. Lett., 2019, 10, 5552–5556.

---