Vibrational Response of Felodipine in the THz Domain: Optical and Neutron Spectroscopy Versus Plane-Wave DFT Modeling

We present a joint experimental and computational terahertz (THz) spectroscopy study of the most stable polymorph (form I) of an antihypertensive pharmaceutical solid, felodipine (FLD). The vibrational response has been analyzed at room temperature by combining optical (THz-TDS, FT-IR, THz-Raman) and neutron (INS) terahertz spectroscopy. With the challenging example of a large and flexible molecular solid, we illustrate the complementarity of the experimental techniques. We show how the results can be understood by employing ab initio modeling and discuss current progress in the field. To this end, we employ plane wave formulation of density functional theory (plane wave DFT) along with harmonic lattice dynamics calculations (HLD) and ab initio molecular dynamics (AIMD) simulations. Based on a comprehensive theoretical analysis, we discover an inconsistency in the commonly accepted structural model, which can be linked to a distinct librational dynamics of the side ester chains. As a result, only a moderate agreement with the experimental spectra can be achieved. We, therefore, propose an alternative structural model, effectively accounting for the influence of the large-amplitude librations and allowing for a comprehensive analysis of the vibrational resonances up to 4.5 THz. In that way, we illustrate the applicability of the computationally supported THz spectroscopy to detect subtle structural issues in molecular solids. While the provided structural model can be treated as a guess, the problem calls for further revision by means of high-resolution crystallography. The problem also draws a need of extending the THz experiments toward low-temperature conditions and single-crystal samples. On the other hand, the studied system emerges as a challenge for the DFT modeling, being extremely sensitive to the level of the theory used and the resulting description of the intermolecular forces. FLD form I can be, hence, considered as a testbed for the use of more sophisticated theoretical approaches, particularly relying on an advanced treatment of the van der Walls forces and going beyond zero-temperature conditions and harmonic approximation.

A majority of active pharmaceutical ingredients (APIs) are distributed in solid forms, where the extent and the type of the molecular ordering and the associated intermolecular forces stay at the source of variation in their physicochemical properties, particularly the dissolution rates, which determine their bioavailability. This can be associated with the molecular forces, binding the molecules together and restricting their mobility. Since the phonon modes manifesting in the THz spectra mainly originate from external, lattice vibrations, they are very sensitive to crystal packing and so the intermolecular interactions, facilitating a better understanding of the crystalline compounds at the molecular level. This makes THz spectroscopy one the most sensitive techniques for the identification and characterization of polymorphic pharmaceutical solids.
Exploration of the terahertz range of frequencies can be realized utilizing absorption and scattering vibrational spectroscopy techniques. The first terahertz absorption spectra were reported a couple of decades ago by using Fourier-transform infrared spectroscopy (FT-IR). A completely new strategy to study molecular excitations in the far-infrared regime has been established with the development of the time-domain spectroscopy (TDS), which currently became the standard method of terahertz spectroscopy [22,23]. FT-IR has been superseded by THz-TDS mainly because of its speed, brighter sources, and higher signal to noise ratios, thanks to its direct access to complex dielectric properties. These make THz-TDS more appropriate for the transmission measurements. Complementing the absorption spectroscopy, a range of scattering techniques, including the inelastic scattering of light (Raman), X-Ray (IXS), and neutron (INS), can be used to probe the vibrational dynamics in the THz domain.
Tremendous interest in terahertz spectroscopy entailed the need for a robust computational strategy, which enables understanding of the complex information provided by the THz experiments. The complexity of the low-frequency vibrations and their sensitivity to the structural environment call for the use of advanced theoretical approaches. This can be satisfied by employing ab initio calculations, where a leading first-principles approach that can yield reliably simulated spectra is density functional theory (DFT). Notwithstanding the progress made on the simulation of the solid-state THz response, one should underline that this remains one of the greatest challenges for computational molecular spectroscopy, where no universal approach exists to date. The description of the THz spectra of crystalline compounds requires taking the periodic arrangement of atoms into account, which, however, can be realized in many ways with different solid-state DFT codes. However, each code implements the formalism in a different way, putting in doubt the reproducibility of such predictions. The problem of reproducibility in DFT calculations of solids has been recently extensively discussed by Lejaeghere et al. [24]. Irrespective to selected computational strategy, the phonon calculations in the THz regime give rise to many issues with convergence, and a systematic step-by-step approach is necessary to achieve well-converged results. The energy expansion in phonon calculations stems from the principal assumption that the system is in mechanical equilibrium and so that all atomic forces are zero. Exceptionally high accuracy in the description of the Hellman-Feynman forces is an essential condition to accurately explore the shallow-and low-energy parts of the multi-dimensional potential energy surface (PES). This makes the calculations computationally expensive.
In this work, we briefly discuss the trends and advances in the modeling of the terahertz spectra and present the case study of felodipine (hereafter, FLD). FLD (Renedil, Plendil) is a popular, second-generation dihydropyridine calcium channel blocker, commonly used as an antihypertensive agent. With the research example of a large and complex pharmaceutical solid, we illustrate how modern ab initio solid-state modeling can support the complementary experimental study. To this end, we employ a combined experimental protocol, involving optical (THz-TDS, FT-IR, THz-Raman) and neutron (INS) terahertz spectroscopy. In our previous work Pajzderska et al. [25], we have extensively characterized the structural properties and molecular dynamics of the title system, with an emphasis on the use of computationally supported neutron scattering and nuclear magnetic resonance (NMR) techniques. Here, the joint experimental and theoretical study is presented aimed at the interpretation of the THz vibrational resonances observed by each technique.
The paper is organized as follows. We first discuss the experimental and theoretical principles of optical and neutron terahertz spectroscopy and present the computational strategies used to facilitate interpretation of the spectroscopic measurements. We next present the case study. While starting from a general discussion of the terahertz spectra, we turn into the analysis of the results with the help of ab initio modeling, providing a comprehensive interpretation of the experimental results. Following-up our previous findings [25], we discover an inconsistency in an idealized model of the crystal structure [26]. Based on the combined experimental and theoretical protocol, we build an alternative, guess model of the crystal structure, allowing for an improved description of the experimental THz data.

Computational Details
All theoretical calculations presented here were performed in periodic boundary conditions (PBC) with the plane wave/pseudopotential formulation of DFT implemented in CASTEP code (version 16.1) [27,28].
We based on the model of the polymorph I (monoclinic, space group P2 1 /c), originally solved from a single-crystal diffraction data set by Fossheim (CCDC: DONTIJ) [26]. In the present work, we rely on the fixed-cell calculations, where the cell parameters were derived from the Rietveld refinement to the room temperature powder X-ray diffraction (PXRD) data [25]. While finding the resulting spectra inconsistent with the experimental data, we further analyzed some alternative structural models, with the molecules packed within the above-defined unit cell. This has been discussed in details in the main text.
Exchange and correlation were approximated with a generalized gradient approximation (GGA)-type functional from Perdew-Burke-Ernzerhof [29]. We have earlier proven such a strategy to be robust for the prediction of the terahertz spectra for complex molecular systems [30][31][32][33], including similar drug, Lacidipine [34]. Furthermore, the use of the room temperature cell volume was dictated by an occurrence of some discrepancies in phonon predictions when using the fully relaxed cell parameters [25], indicating a non-negligible influence of the temperature and calling for accounting for the cell expansion. In addition, a number of different approximations to the exchange-correlation functional were used, including pairwise, semi-empirical dispersion corrections (DFT-D) [35] from Tkatchenko and Scheffler, vdW(TS) [36].
In all cases, the core electrons were described by the set of hard, norm-conserving pseudopotentials, while the electronic wave functions were defined using a PW basis set, with a kinetic energy cutoff of 1050 eV. The Monkhorst-Pack grid was maintaining the constant k-spacing of 0.05Å −1 . Prior to vibrational analysis, the crystal structures were accurately relaxed to minimize the residual atomic forces. The convergence criteria, adopted all throughout this work, in the variation of the total energy, maximum force, external stress (in the case of full-cell relaxation), maximum displacement and the self-consistent field (SCF) were defined as 5 × 10 −10 eV/atom, 1 × 10 −5 eV/Å, 0.0001 GPa, 1 × 10 −6Å , and 1 × 10 −12 eV/atom, respectively.
The harmonic lattice dynamics (HLD) calculations were conducted in the framework of density functional perturbation theory (DFPT) [37][38][39]. In that approach, the phonon frequencies, ω mq , and the phonon eigenvectors, U mq (κ β), at a given q-vector in the Brillouin zone can be obtained by solving the dynamical equation: [40] where κ denotes all the atoms in the unit cell, α and β stand for the cartesian coordinates, andD q α,β (κ, κ ) is the dynamical matrix, which is the fourier transform of the mass-weighted interatomic force constant matrix (IFC),C q α,β (κ, κ ). The IFC matrix is related to the second derivatives of the energy with respect to the ionic displacements from the equilibrium positions (the Hessian matrix).
The provided phonon frequencies and eigenvectors were used directly for the INS modeling. The calculations of the dynamic structure factor, S(Q, ω), were performed with the help of the AbINS program [41].
The IR/THz activities were modeled also from DFPT, as relying on the electric field perturbations. These required the knowledge of the Born effective charge (BEC) tensor and the dielectric permittivity tensors, where both can be expressed as the 2nd order energy derivatives with respect to external electric field. For a given atom, κ, BEC can be linked either to the change in polarization, P, along a direction, α, according to the periodic ionic displacement, τ κ,α , or to force induced on it in the α direction by the macroscopic electric field, ε β [40].
where 0 is the volume of the unit cell. The mode oscillator strength tensor, S m,αβ , defining the IR/THz absorptivity of a given phonon at the Gamma point (q=0), can be written as: The frequency-dependent dielectric function can be then calculated using the earlier obtained phonon frequencies and eigenvectors and the BEC tensor as the following [40]: where ∞ αβ is the high-frequency dielectric constant (electronic polarization). ∞ αβ is calculated using polarization induced in the α direction when a field is applied in the opposite direction, β, incorporating the effect of electron screening within a solid.
The above-given statements assume that transverse optical modes (TO) are considered. However, in the long-wavelength limit (q − → 0), the phonons breaking the center of symmetry are associated with a macroscopic polarization, leading to a longrange dipole coupling. The macroscopic electric field generated in that way exerts a force on the atoms, thus affecting the phonon frequency. To account for this, an effective dynamical matrix needs to be build by adding a non-analytical contribution, NAC q α,β (κ, κ ), to the IFC matrix. This contribution links the BEC and the ∞ αβ tensors with the applied electric field direction.
In this work, the related effects arising from the macroscopic electric field were included in the THz-TDS simulations with the help of the PDiel code [42]. PDiel uses an effective medium theory to account for an influence of the induced electric field on the complex frequency-dependent permittivity tensor of the crystal and furthermore examines the influence of the supporting medium and the sample morphology on the resulting absorption spectrum.
The evaluation of the Stokes Raman intensities requires a third-rank tensor. This was done numerically for the TO modes with the use of a hybrid DFPT/finitedisplacement scheme, by differentiating the electric polarizability, χ = ( ∞ αβ −1)/4π , with respect to the displacement coordinate: where C is the prefactor transforming predicted activities into the intensities to be directly compared with the experiment. The prefactor C accounts both for the laser excitation wavelength, ω laser , as well as the temperature, T: where the Bose factor, n(ω i ), is further defined as In addition to the HLD calculations, the classically thermostatted (NVT; Nosé-Hoover), time-evolved simulations were performed in the framework of Born-Oppenheimer molecular dynamics (BOMD) scheme, in order to account for the finite-temperature effects. The numerical conditions were kept the same apart from the SCF convergence accuracy, which was reduced down to 1 × 10 −7 eV/atom. After successful equilibration, the production run of 12.5 ps was collected. Using 320 processors at the PL high-performance computing facility PROMETHEUS, [43] these calculations took about 5 days per ps.
For the MD simulations, the collected trajectory served as an input for modeling of the hydrogen-projected vibrational densities of states (VDOS) through the calculations of the velocity autocorrelation function (VACF), allowing for a direct comparison with the INS experiment. The GDOS was calculated by Fouriertransform of the velocity autocorrelation function (VACF) provided by the AIMD simulations: where ν(t) κ is the velocity of the species κ at time t.
The collected trajectory was also used for calculations of the anharmonic THz-IR spectrum. To this end, the dipole moments, μ, under periodic boundary conditions were approximated at each step of the reference trajectory by means of the maximally localized Wannier function scheme [44][45][46][47]. The THz-IR spectrum was obtained from the VACF of the dipole moments approximated in that way. The localization of Wannier functions was performed with CP2K [48,49], treating the CASTEP trajectory as an external input. PBE in combination with plane wave (500 Ry) / gaussian atom-centered double-zeta quality basis sets (DZVP-MOLOPT) was used along with the Goedecker-Teter-Hutter (GTH) pseudopotentials. The analysis and the post-processing of the production runs were made with the help of TRAVIS code [50,51].

Vibrational Spectroscopy
In the present study, all experiments were performed for a hydrogenous sample of FLD (Cipla Ltd., Mumbai, India). The sample has been previously extensively characterized through calorimetry, PXRD, 13 C CP/MAS NMR, 1 H NMR relaxometry, and neutron scattering techniques [25].
Standard THz time-domain spectroscopy (THz-TDS) has been used for the determination of the dielectric function. Teraview TPS 3000 unit with accessories in transmission configuration was used. The main parameters of the system are an accessible spectral range of 0.06-3.5 THz, signal to noise greater than 4000:1, and the spectral resolution better than 0.06 THz. The chamber was purged with dry air to eliminate water vapor. FLD was mixed with a high-density polyethylene (HDPE) powder to reduce the attenuation from the sample. The material was ground using a mortar and pestle to reduce the particle size and to avoid the scattering loss. Next, the sample was mixed with the PE matrix powder. The sample was then pressed into a pellet (400 mg weight and 13 mm in diameter), using a hydraulic press. As a reference, a pellet made of pure polyethylene, with the same diameter and the weight of 360 mg, was prepared. The spectrum was acquired by the accumulation of 3600 scans and transformed into the absorption form using a standard, 3-term Blackman-Harris apodization function.
In order to extend the spectral range toward higher-frequencies, we have employed Fourier-transform infrared spectroscopy (FT-IR). The measurements were performed using a dry air-purged Thermo Scientific Nicolet iS50 spectrometer, equipped with a DLaTGS detector and a Globar IR-source. The measurements were performed in an attenuated total reflection mode (ATR; using a monolithic diamond crystal). The spectra were recorded with a resolution of 0.06 THz, by accumulating 32 scans.
The Raman scattering measurements were carried out in the backscattering geometry, using a Renishaw inVia spectrometer, equipped with a confocal DM 2500 Leica optical microscope, a thermoelectrically (TE) cooled CCD detector and an Ar + ion laser, working at 514.5 nm. In the standard setup, the edge filters were used (Rayleigh cutoff at 3 THz). The spectral range was extended down to 0.3 THz by using a NExT filter from Renishaw (composed of a double monochromator, working in Czerny-Turner configuration). In both cases, the instrumental resolution was better than 0.06 THz. Twenty accumulations were collected, where the exposure time was set to 20 s at the laser power of 1 mW. The position of Raman peaks was calibrated before collecting the data using a Si sample as an internal standard.
The inelastic neutron scattering (INS) experiment was performed using a coldneutron time-of-flight spectrometer IN6, located at the European neutron source Institut Laue Langevin (ILL) in Grenoble (France) and presented elsewhere [25]. The spectrum was recorded in the energy-gain mode, with the incident wavelength of 5Å. The data were processed to obtain the dynamic structure factor S(Q, ω) that was finally converted into the generalized phonon density of states (GDOS; G(ω)), with a Q-cutoff of 1.5Å −1 , following standard conventional mathematical relations: where B(ω, T ) is the Bose factor. GDOS accounts for the temperature dependence of the S(Q, ω), related to the thermal population of phonon modes.

Current Trends and Advances in Modelling of Terahertz Spectra
The accurate calculations of the optical THz spectra are known in the literature since the last decade. This has been stimulated by the development and implementation of the perturbative schemes, accounting for the interactions with an external electric field [37][38][39][52][53][54][55].
In principle, the major difference in the implementation of DFT formalism in solid-state codes arises from the selection of the definition of the basis for the single-particle wave functions, used to solve the Kohn-Sham equations and related numerical issues. In solid state, the wave function can be expressed by a periodic Bloch function expanded in terms of some set of mathematical basis functions and multiplied by a complex phase factor, whose wavevector is drawn from the first Brillouin zone of the reciprocal lattice of the crystal [56]. Several different approaches have been developed for the use in solid-state codes, including plane waves (e.g., CASTEP [27,28], QUANTUM ESPRESSO [57,58], ABINIT [59]), augmented plane waves (e.g., VASP [60], GPAW [61]), and localized numerical (DMOL 3 [62], FHI-AIMS [63], SIESTA [64]) and Gaussian (CRYSTAL [65]) basis sets or a hybrid gaussian plane wave scheme (CP2K [48,49]) to name just a few.
There are many literature examples of successful use of each of the aforementioned codes, where CASTEP, CRYSTAL, and VASP are perhaps the most popular in the THz community. Historically, the pioneering works by Allis and Korter were devoted to THz-TDS spectroscopy of explosives and drugs. In these papers, the IR activities were approximated numerically through the use of the so-called differential dipole method (using DMOL 3 ) [66][67][68][69][70][71][72]. Soon after, this evolved into the analytical calculations based on the coupled-perturbed Kohn-Sham (CPKS) scheme, using CRYSTAL code [73][74][75], which is continued to date by Korter and co-workers as well as by other research groups [5-7, 76-87, 87-99]. A similar strategy, with the use of the localized orbitals, was presented later the Schmuttenmaer's group, using SIESTA [100][101][102] and Gaussian codes, respectively [103].
An intrinsically different strategy was presented in 2007 by Jepsen and Clark with the seminal paper on the use of plane wave density functional perturbation theory (PW-DFPT) for analytical calculations of optical terahertz spectra [104]. While introducing the use of the analytical response scheme by Gonze and Baroni [37][38][39][52][53][54][55], this paper provides a proven numerical recipe for an accurate modeling of low-wavenumber phonon frequencies. This strategy has been followed with a number of examples presented in the literature, using CASTEP [30][31][32][33][105][106][107][108][109], QUAN-TUM ESPRESSO [106,[110][111][112][113][114], ABINIT [115][116][117], or VASP [118,119], and has been used in the present study. Instead of advocating any of the presented strategies, we further highlight the main differences and limitations of each approach as both have their pros and cons. One should note that among many of the solid-state codes available, the use of plane wave basis sets is the most widespread [56]. In contrast to the localized basis sets (centered at the atomic positions), plane wave basis sets (i) are spatially unbiased; cover the whole cell volume, being free from basis set superposition errors (BSSE); and finally (iii) allow for a systematic convergence, since their size is controlled by a single parameter, the cutoff energy E cut . However, as recently shown, it does not ensure reaching the basis set limit [120]. The main disadvantages in comparison with the atom-centered bases include [56] (i) large numbers of required basis functions, which easily make the calculations inefficient; (ii) serious problem for the evaluation of the Fock exchange term as required by hybrid exchange-correlation (XC) functionals, and (iii) difficulty to describe the electronic structure near the nuclei and therefore the need of using pseudopotentials, which definition can be critical in all the aspects of the computations performed, including lattice dynamics.
While the atom-centered basis sets are superior regarding these aspects, they, however, suffer from a number of different problems, most importantly, from the BSSE issues due to an incompleteness of the basis set as well as from the so-called pseudo linear-dependence problem. In contrast to plane waves, the size and the quality of the basis set do not depend on a single parameter; therefore, its optimization is highly non-trivial. Usually, the high-quality quantum chemistry basis sets in molecular modeling are augmented with diffuse functions, i.e., functions that have smaller orbital exponents than those normally used. These are particularly useful in the cases where the charge distribution is expected to be considerably more diffuse than in the neutral atom, e.g., for the negatively charges species or polar systems, carrying an excessive negative charge in some molecular parts. However, their use in solid-state modeling is generally prohibited. Diffuse components (i) dramatically increase the number of integrals to be explicitly calculated, (ii) call for ultra-high numerical accuracy to avoid the linear-dependence issues, and finally (iii) may lead to overlapping with the functions of other atoms in densely packed solids.
Another issue limiting formulation of universal, gold-standard computational methodology for the modeling of terahertz spectra stays at the heart of density functional theory. Since the precise functional dependence of the energy on the density is not known, DFT relies on the crucial approximation. While generally, the whole functional ZOO is available in the molecular modeling software, very few approaches are accessible in the solid-state codes. The aforementioned problems in calculations of the Fock exchange term with plane wave DFT practically limit itself to the use of the semilocal approximations. On the contrary, the exact Hartree-Fock exchange can be naturally and efficiently introduced via the localized basis sets and this has been exploited in several codes, including CRYSTAL [65], CP2K [48], and FHI-AIMS [63], with a great number of more advanced approximations to the XC correlation functional available. However, as concluding from recent benchmark studies testing a broad range of functionals in solids (LDA, GGA, meta-GGA, and hybrids), overall, none of the approximations can be considered as a universal method that is better than the others when applied to a broad range of solids [121]. The accuracy of XC approximations in the prediction of the harmonic frequencies in the THz domain also remains an open question. It should be underlined here that an excellent agreement between the experimental and harmonic frequencies obtained using DFT calculations is often the result of a fortuitous cancellation of errors. The calculated harmonic frequencies should not agree with the experiment, because anharmonicity can produce shifts in the fundamental frequencies. This has been discussed in details, e.g., by Gill, proposing the approach (EDF2) explicitly designed to yield accurate harmonic frequencies as benchmarked against accurate CCSD(T) golden standard [122]. It however appears that the anharmonicities in the crystals below 200 cm −1 are relatively small [123].
One of the most significant shortcomings in DFT is an inadequate treatment of weak van der Waals (vdW) interactions, stemming from the long-range electron correlation effects. Even in the short to medium range, most semilocal density functionals fail in describing the dispersion forces. A common approximation used for weakly bound molecular crystals (also adopted in this work) is to assume that the most important error caused by omitting long-ranged dispersion is the cohesive stress term and to optimize and compute phonons with the lattice vectors according to diffraction findings. This is done under a fictitious external pressure, mimicking the conditions arising from the cumulative sum of the long-range dispersion interactions. However, great progress was made in the last two decades to tackle this problem by developing a number of schemes for dispersion corrections to the XC functional (DFT+D). This has been discussed in several review papers [124][125][126][127][128][129]. In analogy to the famous "Jacob's ladder" concept for the XC functionals [130], Klimeš used the term "stairway to heaven" to classify and group DFT-based dispersion correction schemes. At the ground level are the methods which do not describe the long-range asymptotics. Simple pairwise correction schemes (e.g., DFT-D2 [131]) sit on the first step; on the second step are approaches that utilize environment-dependent terms (e.g., DFT-D3 [132], vdW(TS) [36]). The long-range density functionals (e.g., vdW-DF [133]) sit on the third step, and at the top are the approaches which go beyond the pairwise additive determinations of dispersion (e.g., MBD [134,135]) [125].
Advanced dispersion corrections greatly improve the description of the structural and cohesive properties of molecular crystals [136]. Extensive benchmarks have proven that the many-body dispersion schemes are particularly promising, allowing for reliable predictions of the polymorphic energy landscapes close to the chemical accuracy [137]. Currently, most of leading solid-state codes give access to the highest steps of the stairway, with, however, limited functionality. Due to increasing complexity, these were found to be particularly difficult to implement for calculating the atomic response properties via the perturbative schemes, limiting their use for modeling of the THz spectra. The problem has been recently solved in ABINIT [138], fully implementing the PBE-D3 schemes, which are also the most advanced ones available in CRYSTAL code, becoming a method of choice in the THz community. CASTEP allows using pairwise dispersion corrections (e.g., DFT-D2, vdW-TS) to this end. At present, however, none of the solid-state codes allows for the use of the many-body dispersion (MBD) scheme for the prediction of the atomic response with the perturbative scheme, whereas this approach was reported as very promising for modeling in the THz domain [139]. In this regard, the phonon frequencies and eigenvectors can be obtained using standard, finite-displacement approach to lattice dynamics. This information is generally sufficient to model the inelastic neutron scattering data, providing a unique way of testing the performance of more advanced theoretical approaches [140]. Recently, several authors remarked on a very good performance of the VdW-DF family of functionals of Langreth et al. [133] in the modeling of the INS spectra [140,141]. Similar remarks can be derived from the results of the THz predictions from the Schmuttenmaer's group [100][101][102].
An alternative approach to lattice dynamics calculations is to abandon the concept of a normal mode and explore the multi-dimensional PES utilizing the molecular dynamics simulations. Such a strategy allows one to go beyond the zero temperature and naturally account for mechanical anharmonicities and thermal effects. This might be of particular importance, since a vast majority of the experiments are performed at room temperature, where the molecular mobility and relaxation effects cannot be ignored. MD allows accounting for mode couplings, which is particularly important in the case of complex, floppy molecules. Furthermore, it gives a natural opportunity to account for molecular disorder in solids at both local and global scale. Applicability of Born-Oppenheimer molecular dynamics (BOMD) simulations in the modeling of the THz spectra was discussed, e.g., in [123]. Recently, Ruggiero and Zeitler have reported successful attempts to combine CRYSTAL and CP2K to account for the anharmonicity and elucidate the nature of phase transitions in several solids [87,91,142]. In this work, we present a similar strategy, combining CASTEP and CP2K to predict INS and THz-TDs/IR spectra of the title system.

The Structural Model of Felodipine
FLD is a long-acting 1,4-dihydropyridine calcium channel antagonist, marketed in a solid form. The molecular geometry of FLD is displayed in Fig. 1a. There are four polymorphs of FLD known to date [143,144], where the monoclinic form I (P2 1 /c; Z = 4) is the most stable one (see Fig. 1b) [26]. While the single-crystal structure was solved at 150 K [26], we rely here on the Rietveld refinement to the room temperature powder X-ray diffraction (PXRD) data for the studied sample [25], resulting in the following cell constants: a = 12.18Å, b = 12.24Å, c = 13.52Å, and β = 116.05 o In all the solid forms, the FLD molecules adopt a similar conformation, where two planar rings (dichlorophenyl, PhCl 2 , and 1,4-dihydropyridine, DHP) are perpendicular to each other. These constitute a rigid molecular framework [25,143,144], whereas the side chains were found to be more sensitive to molecular packing.
In our previous work, we have quantitatively analyzed the intermolecular forces driving its crystal packing [25]. It was found that the stacking forces between the rings, binding the rigid molecular frameworks together, is the most prominent intermolecular interactions, dominating over a medium-strength NH· · ·O bonding. Since the THz spectroscopy techniques primarily explore the external lattice vibrations, The projection of the unit cell of FLD form I (P2 1 /c; Z = 4) according to the refined X-ray diffraction data [25] they allow probing the motions of the molecules confined in the potential of these intermolecular forces.
While no thermal factors were provided by Fossheim [26], we also found that the FLD molecules manifest a prominent thermal-induced dynamics, where the internal molecular mobility has been primarily ascribed to the flexibility of the side ester chains and the methyl groups, which was extensively analyzed with the help of 1 H NMR relaxation and quasi-elastic neutron scattering (QENS) [25].
The analysis of the INS spectra up to ca. 10 THz suggested that the floppy molecular fragments are described by a considerably anharmonic vibrational behavior induced by large changes in the amplitudes of atomic motions. While in this case the INS spectroscopy is primarily sensitive to hydrogen motions, it does not tell much on the vibrational response of other atomic species. Therefore, we extend the analysis toward the non-hydrogenous species by employing the THz-TDS/IR and Raman, which are primarily sensitive to displacements of heavier atoms. In that way, we intend to provide complementary information on the THz response of the title drug by referring to both flexible side fragments as well as the rigid molecular framework.

Experimental Spectra: General Considerations
In this work, we employ an experimental protocol which combines INS and optical THz techniques (THz-TDS, IR, Raman). The resulting spectra are presented in Fig. 2.
Large molecular unit of FLD gives rise to very complex THz spectra due to a great number of vibrations involved. In total, 528 fundamentals are to be observed in the vibrational spectra of the title system. According to the single-crystal X-ray diffraction (SXRD) data [26], the crystal structure is centrosymmetric, and so the mutual exclusion rule prohibits the modes to be both infrared and Raman active. In brief, Fig. 2 The room temperature terahertz spectra of FLD form I presented in the range of 0.25-4.5 THz. The THz-TDS (red) and THz-Raman (blue) spectra have been combined with the ATR-FT-FIR and normal Raman spectra, respectively (merged at 3.75 THz). The room temperature INS spectrum (black bold line) is compared to the ones measured at 250, 200, and 150 K (the lines narrowing upon cooling), respectively, to expose the smeared bands analyzed. The peak fitting was performed with the set of Lorentzian functions the reducible representation of the monoclinic form I of FLD can be decomposed as follows: where 1A 1u and 2B 2u are the acoustic modes with a zero frequency at the point (q = 0). All gerade modes are Raman active, whereas the ungerade ones are THz/IR active. Due to a low symmetry of the system, no silent modes are present. According to the -point phonon calculations (fixed-cell PBE), the discussed spectral range up to 4.5 THz covers 71 non-zero-frequency modes. These were analyzed in details in Fig. 3, presenting the mode classification according to their internal, rotational, and translational character. The internal vibrations probed in this range involve complex torsional deformations, considerably delocalized over the internal coordinates of the whole molecule. An inspection of the figure shows that the mixed, internal/external nature of vibrations is manifested all throughout the frequencies covered by the THz-TDS and THz-Raman experiments. The external modes become more prominent below 2.5 THz.
As stemming from the dipole approximation, the THz-TDS/IR and Raman techniques are subject to optical selection rules. On the contrary, INS spectroscopy relies on the interaction with the nuclei not the electron cloud. Therefore, it provides a  [26]. Stacked histograms showing energy decomposition arising from the contributions from the translations and rotations of the molecular center of mass (CM) and the internal molecular vibrations (see [106] for a detailed description of the numerical procedure) microscopic insight into a material in the absence of selection rules, where, in principle, all modes are observable [145]. The intensity of an INS band is defined by the so-called scattering law, S(Q, ω). For a powdered sample, the spectral intensity of a phonon with frequency, ω, is proportional to a simplified expression: where n is the quantum number of an ith mode (n = 1 for a fundamental; n > 1, 2, 3, and so forth for overtones and combination bands) and σ (barn) stands for the scattering cross section. Q (Å -1 ) denotes the momentum transfer. U i is the mean square atomic displacement in the considered vibration, and the exponential term is known as Debye-Waller factor (DWF). U T ot stands for the total root mean square displacement of all the atoms in all (both internal and external) modes, and its magnitude is partially determined by the thermal motion of the molecule. DWF results in damping of the observed intensity. This can be reduced by cooling the sample and so the spectra are typically recorded at the cryogenic conditions [146]. The Debye-Waller term also depends on the momentum transfer, so it can be reduced by recording the spectra at low Q values (as presented here), similar to standard IR and Raman spectroscopy which exclusively probe the point. A review of the use of INS in molecular spectroscopy can be found, e.g., in [147], whereas some recent advances were discussed in [148]. Depending on the way how the energy of the scattered neutrons is analyzed, the INS experiments can be realized either with the use of the so-called direct-(present case) or indirect geometry spectrometers. Advances of using each type of instrument are discussed in details in [146]. Since Eq. 11 is purely mechanical, it makes the computational predictions straightforward as the product of the neutron scattering cross section and the amplitude of vibration. On the contrary, IR and Raman techniques rely on the interactions with an electron cloud, making the response to nuclear dynamics indirect. Moreover, in Raman spectroscopy, the polarizability for each atom normally grows with the number of electrons. Therefore, it is rather insensitive to hydrogen dynamics. While IR spectroscopy is sensitive to electric dipole, the activities are dominated by a considerably charged species. This is roughly reflected in Fig. 4, giving the percentage ionic contributions to the total atomic displacement in each mode under interest. The strongest contributions from hydrogen reflect in the most intense bands observed experimentally at ca. 2 and 3.5 THz. The visible contributions from oxygen spread all throughout the studied range, express both external motions and the ester chain deformations, whereas the contributions from the chlorine atoms follow the distribution of the Raman intensities observed experimentally and can be attributed to both external modes and the intramolecular PhCl 2 deformations.

Theoretical Spectra: Insights from the Crystallographic Model
Among the 71 predicted modes, 34 are allowed to be observed with THz-TDS/IR and 37 ones are Raman active. This makes the theoretical description challenging. The theoretical IR spectra are displayed in Fig. 5. The theoretical predictions provide only a qualitative match to the experimental data, and the assignment is rather nonobvious. The source of the mismatch can be ascribed either to (i) the shortcomings of the theoretical methodology applied or (ii) to the use of an improper structural model. Here, we consider the former issue, whereas the structural aspects are discussed in the next section. Figure 5a examines the influence of the use of alternative generalized gradient approximations. Most GGA functionals can be classified in terms of an analytic function, known as the exchange energy enhancement factor, F x (s), used to improve over local density approximation (LDA). One may consider the PBE approximation as the standard one as obeying many of the exact theoretical conditions and devoiding any  [26]. Stacked histograms show partial contributions to the total atomic displacement for each phonon mode  [26] fitting parameters. The so-called "soft" revision, dedicated to solids, PBEsol [149], introduces F x (s) that increases more slowly with the reduced density gradients. On the contrary, we also consider a "hard" GGA formulation, rPBE [150]. Finally, the PBE augmented with the pairwise van der Walls corrections from Tkatchenko and Scheffler is denoted as PBE-TS. The presented results indicate a dramatic influence of the selected GGA approximation. While these results were obtained at the cell vectors fixed to the room temperature data, the influence of the cell relaxation was examined with the help of the PBE-TS (see the FOPT). Whereas none of the approaches was found to be superior, we further limit ourselves to the use of the most standard one (PBE). In this regard, the system should be considered as a candidate for the future benchmarks using more sophisticated ab initio approaches, particularly relying on more advanced corrections for the dispersive interactions combined with hybrid functionals. Figure 5b presents the results of the fixed-cell, HLD PBE predictions as comparing the spectrum generated from the traverse-optical modes (TO) with the ones taking into account the supporting medium and the particle shape using an effective medium method. [42] This was done by considering the spherical and plate particles embedded in the HDPE matrix. The plate shape was selected since such a morphology was observed both in the single-crystal experiment [26] and with a microscope (micro-Raman experiment presented here). The three families of the [h k l] faces are considered here as giving the most distinctive spectra among all the combinations tested. Nevertheless, from the inspection of the figure, it appears clear that the molecule is not much polar and the long-range dipole coupling is not particularly important in this case. Therefore, we do not link the discrepancy to the influence of the sample morphology.
In Fig. 6, we finally compare the results of standard PBE calculations with the full collection of the experimental results. However, since only moderate agreement between both sets of data has been achieved, it does not allow for an unambiguous interpretation of the recorded spectra.
The harmonic lattice dynamics predictions presented in the figure have been further compared to the results of ab initio MD simulations. Here, we limit the discussion to the IR and INS data, since the MD simulations of the polarizability evolution was beyond our computational capabilities. While comparing both sets of data, one may note some non-obvious frequency shifts when accounting for the temperature effects, whereas most of the predicted bands are blueshifted with the frequency. At this stage, it is, however, questionable to which extent this trend is physical, reflecting the influence of the temperature and dynamics, or if it is biased by an increasing artificial internal stress in the fixed-cell NVT simulations. The MD strategy was dictated by omitting of the dispersion corrections, whereas ideally, the NVT runs should follow NPT equilibration, which makes no sense without the use of the vdW-corrected schemes. By inspection of the figure, one may, however, note that the inclusion of the temperature effects improves the agreement with the experimental spectra in terms of the intensity distribution, pointing to the importance of the related dynamic factors. This is further suggested by a broadened nature of the experimental bands, suggesting that the vibrational relaxation and the mode coupling effects may be of importance in this case.
One may question, why the above described MD simulations do not properly tackle the dynamic contributions from the ester group librations. One of the causes might be a short time scale of the performed simulations, as well as the limited box size. However, more importantly, the Born-Oppenheimer MD does not account for the zero-point energy (ZPE) contributions. This can be overcome by employing quantum thermostatting, which its use in the simulations of the THz regime is, however, limited [151].
While only a qualitative description of the experimental spectra is possible, the analysis of the mode eigenvectors shows that most of the IR activities arise from the activation of the ester group deformations, the dichlorophenyl ring torsions and waggings, and the motions affecting the hydrogen bond, which by its nature leads to pronounced dipole fluctuations.
While referring to the analysis of the intermolecular interactions [25], we note that the most stabilizing forces in the crystal are due to anti-parallel stacking interactions, binding the rigid ring frameworks together. The lowest frequency modes (below ca. 1.75) can be roughly ascribed to the shearing translation modes, mixed with torsional deformations. Fig. 6 Experimental Raman, THz-TDS/IR, and INS spectra measured at 300K compared to the results of the theoretical predictions by means of harmonic lattice dynamics (HLD at 0 K) and ab initio molecular dynamics (MD) simulations (NVT at 300 K). See the text for more numerical details. All the calculations rely on the crystallographic model proposed by Fossheim [CCDC: DONTIJ] [26]. The arrows and dotted lines, linking the experimental data with the theory, serve as a guide for the eye Furthermore, while analyzing the crystal structure, one can note that a complex network of a medium-strength NH· · ·O hydrogen bonds is formed between the adjacent molecules through the connection with the ester COOCH 3 groups. The hydrogen bonds were found to be non-linear. Therefore, any deformation mode affecting both the N· · ·O distance and the related angle results in strong dipole fluctuations. The most prominent THz-TDS band at 1.87 THz is a by-product of such a displacement. This is corroborated by the upper-frequency shifts revealed by the MD simulations. This most intense THz-TDS signal also involves some contributions from the ring deformations.
The spectral range above 2 THz reflects strong contributions from the ester chain librations, τ [COC 2 H 5 ]. The most intense band predicted theoretically at ca. 3.5 THz does not match the THz-TDS observations. However, the lowering of the intensity was observed through the MD simulations, which point out the importance of the dynamic effects.
Similarly, the MD simulations greatly improve the description of the INS intensity distribution, which is purely due to large-hydrogen displacements, particularly coming from strong librations of the CH 3 and COC 2 H 5 species.

Toward an Improvement of the Structural Model
While discussing the mismatch between the theoretical and experimental spectra, one should consider a possibly incorrect assignment of the polymorphic form. The comparison of the PXRD patterns for all known FLD polymorphs can be found in [143]. We note that the PXRD patterns recorded for the studied sample (see [25]) unequivocally confirm the presence of the form I in our case. While Surov et al. confirmed that the form I is found to be the thermodynamically most stable phase at ambient conditions [143], we further proven that this structure is stable down to 20 K, since no phase transition was detected with the 1 H NMR relaxometry experiments [25]. We also note that FLD is a chiral compound and the form I is a racemate structure. The polymorphism of racemic FLD and the formation of a series of solid solutions in the binary system of its enantiomers were studied in details by Rollinger et al. [152]. An extensive analysis through PXRD and FT-IR in comparison with these data confirms that our sample is a racemate. Taking both factors into account, we conclude that the cell shape, the crystal packing, and the rigid core arrangement adopted in the simulation box are valid, staying in line with the crystallographic model proposed by Fossheim [26].
We, hence, conclude that the main sources of the mismatch between the output of the ab initio calculations and the experimental data can be ascribed to the dynamic effects, modulating the average, idealized crystallographic structure. We, believe that such an inconsistency may be due to an imperfect description of the room temperature structure.
The analysis of the intermolecular interactions and the potential energy surface (PES) indicates that the molecular conformation of FLD is locked by strong intermolecular potential, apart from the COC 2 H 5 fragment, which was shown to undergo strong librations as revealed by the 1 H NMR experiments [25]. Hence, the molecules were found to be dynamically disordered at the OC 2 H 5 sites. While this was not reported by Fossheim at 150 K, we believe that the adopted crystallographic model of the polymorph I is idealized and not fully valid for the description of the room temperature and/or powder samples. Interestingly, most of known FLD polymorphs and solvates solved to date manifest some kind of positional disorder of large atomic displacement parameters at the COOC 2 H 5 sites [143,144,[153][154][155]. This shows that the dynamics of the ester group may be of crucial importance in the characterization of all polymorphs and pseudopolymorphs of FLD known to date.
In order to improve the description of the THz spectra, we have examined a number of alternative configurations, taking the COC 2 H 5 librations into account. We have based on our previous findings from the analysis of the potential energy surface presented elsewhere [25]. With the help of solid-state DFT (PBE-TS), we illustrated the presence of two local minima on the PES landscape, describing the ϕ[COCC] coordinate. These two local minima were described with the energies higher by 9.90 kJ/mol and 25.05 kJ/mol per molecule (PBE-TS), respectively, as referring to the global minimum. Hybrid DFT modeling for an isolated molecule model has revealed that the associated conformational changes are characterized by the energy changes within 4.2 kJ/mol [156]. While the solid-state DFT modeling finds the conformational changes rather unexpected, the excursions toward the metastable configurations were clearly evidenced in the 1 H NMR relaxometry experiments and MD simulations [25].
Following our findings, we have performed the constrained PES scans, with the ϕ[COCC] coordinate varying collectively every 30 o for all four molecules in the cell. This was followed by the relaxation of all remaining internal coordinates (with the associated forces minimized with the accuracy better than 5 × 10 −5 eV/Å) and by subsequent vibrational analysis. In this case, the global cell symmetry was formally turned off (P1), but practically preserved within 0.0001Å, due to the presence of the constrains. Such a procedure allowed us to find the optimized metastable structures (very close to the P2 1 /c point group) and probe the influence of the ester chain orientation on the THz response.
The vibrational analysis revealed that some of the structures can be considered as mechanically stable, where no imaginary phonon frequencies and three, well-defined zero-frequency modes are found. The comparison of the experimental and theoretical optical THz spectra obtained for these models is given in Fig. 7  One may point that the crystallographic and the ϕ[COCC] Eg. +180 o models represent two extreme cases. Possibly, the experimentally recorded spectrum may reflect an intermediate case, where the interconversion between the configurations is very fast and so an unresolved spectrum is recorded. While the redshifts of the most intense bands with respect to the experimental THz-TDS spectrum characterize the former model, the later one exhibits blueshifted distribution of the IR-active modes. Interestingly, an opposite trend can be observed in the case of the Raman-active modes.
The harmonic lattice dynamics calculations allow to analyze the thermodynamic stability of each alternative structure. This has been presented in Fig. 8, where three cases are considered. As expected from the PES analysis [25], the alternative models are much less favorable, with the free energy differences beyond 10 kJ/mol per molecule. While no thermal expansion was taken into account, there is a clear monotropic relation between the studied systems. It leaves us with a puzzling observation and open questions on the extent of the disorder; the role of the configurational entropy; potentially metastable nature of the dynamically disordered structure; and, consequently, the source of the FLD trapping in a clearly less favorable state.
It also puts a question on the sample microstructure and possible formation of domains. This was partially explored by Mishra et al. [157] where an intergrowth polymorphism in single crystals of FLD forms I-III was studied with nanoindentation. While distinct properties for two polymorphic forms within the same crystal of form II were identified, the authors did not reveal any intergrown properties in the polymorphs I and III.  While at this stage we are unable to address these questions, we use the ϕ[COCC] Eg. +180 o as the guess model, allowing to reproduce the experimental spectra to a great extent.

Analysis of the Terahertz Spectra: Insights from the Alternative Model
The conformation of the FLD molecule, extracted from the alternative ϕ[COCC] Eg. +180 o model, is displayed in Fig. 9a. The only visible difference is See the text for more details). Stacked histograms showing energy decomposition arising from the contributions from the translations and rotations of the molecular center of mass (CM) and the internal molecular vibrations (see [106] for a detailed description of the numerical procedure) due to the ester chain orientation. Figure 9b presents a superposition of both models considered. It is clear that the conformational change is accompanied by small displacements of the FLD molecules. The optimized structure can be described by P-1 symmetry (within the tolerance of 0.00001Å) or P2 1 /c (within 0.0001Å). This is reflected by the phonon analysis, where the mutual exclusion rule is nearly perfectly fulfilled. However, the two lowest energy modes are predicted with non-zero intensity for both IR and Raman spectroscopy and have been observed experimentally. This may be an evidence of the temperature-induced symmetry breaking within an imperfect crystal phase, which accompanies the ester group librations. Figure 10 provides the decomposition of the phonon modes into translational, rotational, and internal contributions. Compared to the output of the simulations with the SXRD Fig. 11 Experimental Raman, THz-TDS/IR, and INS spectra measured at 300K compared to the results of the theoretical predictions by means of harmonic lattice dynamics (HLD at 0 K), obtained with fixedcell PBE calculations for the alternative ϕ [COCC] Eg. +180 o model. The arrows and dotted lines, linking the experimental data with the theory, serve as a guide for the eye τ, twisting; δ, bending; ω, wagging; ν, stretching; γ , out-of-plane deformation; lib., libration. DHP and PhCl 2 stand for the rings constituting the molecular core. COOCH 3 and COOC 2 H 5 denote the side ester chains. CH 3 are the methyl fragments model, one can see clear differences, which means that small displacements of the molecular centroids affect the modes to a great extent. In Fig. 11, we finally present the comparison of both sets of the experimental and theoretical spectra. According to the figure, very satisfying match to the experimental spectra has been achieved, allowing for a trustful band assignment. The most prominent experimental bands were then analyzed in details in Table 1. A detailed inspection of the table reveals that the range below 2 THz is highly dominated by the modes of an external character, making their assignment challenging. The analysis of the internal modes contributions gives some insight into the nature of the bands dominating in each spectrum.
In Raman, most of the bands are associated with the PhCl 2 contributions, particularly due to aromatic character of the ring and the presence of electron-rich chlorine atoms. In that way, the technique seems to be very sensitive to the environmental interactions of the aromatic part of the molecular core.
IR response is highly dominated by the torsional motions of the oxygen-rich ester groups. The concomitant C 2 H 3 motions are therefore spread all throughout the THz range. Following the associated frequencies, we see the blueshift compared to the results from the Fossheims's model analyzed in the previous section. This explains the blueshifts in the BOMD while activating the ester chain librations, confirming their physical meaning. As mentioned, any out-of-plane deformation and or external mode resulting in a displacement of the DHP ring affect the hydrogen bonding. The same states for the torsional contributions from the COOCH 3 fragments. Both contribute to the IR intensity by dipole fluctuations.
It appears that in the present case both Raman and INS are hardly sensitive to the hydrogen bond fluctuations. INS is a superb sensitive to the large-amplitude librations of methyl and ethyl groups. The latter contributions make the spectrum more similiar to the THz-TDS one. The fail in the prediction of the frequencies at ca. 3.5 THz, which are clearly associated with the C 2 H 3 librations, suggests that the associated potential is more shallow in reality. Therefore, we may assume that the real structure of FLD form may be in between of the two extremes models considered in this work.

Conclusion
In this work, we emphasize on the joint use of THz-Raman, THz-TDS, and INS spectroscopy to uncover the vibrational fingerprints of large molecular solids. Focusing on the application of plane wave DFT, we extensively discussed its use for the modeling of the terahertz spectra.
We focused on a challenging research example of felodipine form I. The wealth of experimental data collected here facilitates further study of the crystal polymorphism and pseudopolymorphism, as well as supports the analysis of the FLD samples of different origin.
However, this study calls for a future extension toward the low-temperature conditions, by revealing very significant influence of the temperature, where the cell expansion and the dynamic effects are of key importance. We discuss different factors affecting the theoretical predictions. In this regard, we uncover a critical role of the large-amplitude ester group librations, which contribute to the spectra all throughout the THz range.
We consider FLD form I as a testbed, calling for the use of more sophisticated theoretical approaches in the future, dealing both with the van der Walls forces and limitations of the GGA approximation. More importantly, the critical role of the dynamic effects put a quest for the use of more sophisticated MD simulations, accounting for the zero-point energy contributions.
While this is beyond the current state-of-the art of THz modeling, we propose an effective, guess model of the structure, mimicking the influence of the ester chain librations. This allows for a reliable reproduction of the experimental spectra with the help of the harmonic lattice dynamics calculations. In this regard, we provide a comprehensive analysis of the THz spectra. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.