Bulk and Surface Conformations in Solid-State Lovastatin: Spectroscopic and Molecular Dynamics Studies

Abstract

Conformational flexibility in molecules can give rise to a range of functional group terminations at crystal surfaces and dynamic disorder in the bulk. In this work, we explore the conformational behavior of the drug molecule lovastatin in the crystallographically disordered solid and at crystal surfaces through a combination of computational modeling and spectroscopy. Gas-phase and periodic quantum-chemical calculations are used to study the potential energy surface associated with rotatable bonds to examine the disorder in bulk. These calculations are combined with vibrational and X-ray photoelectron spectroscopy measurements to obtain insight into the conformations in bulk and at the surface. Our MD simulations show that the bulk disorder is driven by cooperative motion of the butyl group on the S-butanoate moiety along one crystallographic direction beyond a unit cell. The calculations show that the O-H group can rotate relatively freely between two low-energy conformers in the gas phase but is locked in position by intermolecular H-bonding interactions in the bulk crystal, and we find tentative spectroscopic evidence for the second conformer being present at the surface. We also comment on the relative utility of these different techniques for studying molecular conformation in bulk and at surfaces and highlight possible areas for future developments.

1. Introduction

The conformations of functional groups terminating crystal surfaces can influence important material properties including crystal morphology [1,2], wettability [3,4,5], dissolution rates [6] and adhesive properties [7]. These properties are often particularly important for pharmaceuticals as they affect both the downstream processing and the bioavailability of the end product [8]. Flexible molecules therefore pose a challenge both because they can present functional groups in different orientations at surfaces, and because the presence of multiple conformations in the bulk can influence intra- and intermolecular interactions, crystal packing and give rise to dynamic disorder. A holistic understanding of particle properties requires an understanding of the intrinsic conformational flexibility of the molecule and how this impacts on any bulk dynamic disorder as well as the chemical behavior of functional groups present at the crystal surfaces and the presence of impurities on the true particle surface.

Flexible molecules with multiple bond torsions are capable of forming several “conformers”, defined as molecular structures with distinct topologies that correspond to energetic minima on the structural potential energy surface (PES) [9]. The need to understand the conformational behavior of flexible molecules has led to widespread interest in exploring the energetics associated with conformational changes, understanding the resulting polymorphism [9,10], identifying spectroscopic signatures of different conformers [11,12,13] and applying these insights to crystal structure prediction [14]. Analysis of polymorphs in the Cambridge Structural Database (CSD) indicated that if two topologies result from a difference in a torsion angle Δθ > 120°, they can be considered “conformational polymorphs” [9]. These conformers tend to have energy differences in the order of 10 kcal mol⁻¹, undergo slow inter-conversion in solution and therefore crystallize as distinct polymorphs [15]. If, on the other hand, conformers are separated by a small energy barrier and interconvert rapidly, the conformational flexibility is not likely to affect the polymorphic outcome of a crystallization [15] and may instead lead to the presence of dynamic disorder in a solid.

The structural PES of a molecule in the gas phase can, in most cases, be mapped directly to the conformations observed in the solid state [9] and have been used to identify new conformers and to calculate spectroscopic signatures to help assign experimental data [16,17,18,19,20]. Rotations of, e.g., CH₃, NH₂ and OH groups are usually not considered during conformational analyses, as typically only a thermal average is observed from crystallography. These can, however, play a vital role in the stability of a crystal structure through “conformational adjustments” that facilitate better overall crystal packing in exchange for a small energy penalty.⁹ These changes can be studied using vibrational (e.g., infrared—IR) and nuclear magnetic resonance (NMR) spectroscopy [21,22], and both IR spectroscopy [23] and the NMR nuclear Overhauser effect (NOE) [21,22] have been used to study the conformations of peptides and proteins. Conformational flexibility giving rise to, e.g., molecular reorientation [24,25] and dynamic disorder [26] in organic crystals can be readily investigated using molecular dynamics (MD) simulations, while lattice dynamics modeling can predict the natural thermal motion of atoms and crystal vibrations to provide a point of reference for experimental measurements [27,28]. Such simulations can both help to interpret spectroscopic data [26] and provide enhanced insight into, e.g., the changes to molecular interactions that result from conformational changes due to bond torsions.

In addition to the bulk conformations, there is also the separate question of the conformation of molecules at the crystal surfaces. This is particularly important given that the functional groups exposed at the surface determine the interaction with the environment and therefore control physical properties such as adhesion [7], crystal growth [29] and dissolution rates [30].

We recently demonstrated that polarized Raman spectroscopy can be used to establish the nature and conformation of functional groups at the dominant crystal surfaces of aspirin [27]. Both IR operating in attenuated total reflection (ATR) mode using evanescent waves and Raman using long-wavelength lasers have a penetration depth of microns to mm in an organic crystal, and thus do not provide information about functionality present within nm of the surface. X-ray photoelectron spectroscopy (XPS), which typically probes the uppermost 10 nm of a material and can prove useful when the chemical environment at the surface, is very different to that in the bulk. XPS has, for example, previously been used to establish the conformations of polymers [31] and the protonation states of molecules in the solid state [32].

In this study, we explore the conformational flexibility of the drug molecule lovastatin. Lovastatin (C₂₄H₃₆O₅) is currently used to treat poroketaratosis [33] and hyperchloestrolemia and, more recently, has been explored for the treatment of hypertension. Therapeutically, lovastatin is known to lower the risk of heart disease and strokes in patients with other comorbidities [34] but has very low solubility and bioavailability. An improved understanding of the structure–property relationships in this material is therefore needed to improve its physical properties [35,36] and manufacturability [37,38]. The lovastatin molecule has 65 atoms and eight stereogenic centers and thus can, in principle, exhibit considerable conformational flexibility. Lovastatin crystallizes into the orthorhombic space group P2₁2₁2₁ with four molecules in the unit cell [2,39,40]. Like its analogue simvastatin [26], lovastatin exhibits both positional and dynamic disorder of the S-butanoate arm in the bulk crystal structure [40], which was only recently resolved [2]. In the bulk, we would expect the conformational flexibility to be restricted by intra- and intermolecular interactions. Inelastic neutron scattering studies have identified the presence of stabilizing non-covalent CH$\cdots$HC and CH$\cdots$O interactions in the solid state, but were unable to resolve the dynamics of the O-H bond [34]. On the other hand, molecules at the surface, with fewer restricting interactions, could, in theory, exhibit a broader range of conformations.

In this work, we use a combination of quantum-chemical calculations, MD simulations and computational spectroscopy to explore the conformational flexibility of lovastatin in the gas phase and in bulk. We compare these results to measurements on crystalline powder samples using ATR-FTIR, FT–Raman and near-ambient pressure (NAP) XPS to attempt to shed light on the functional group terminations present on the crystal surface and to examine the conformation adopted by the lovastatin molecules at the true surface and the sub-surface. We also assess the possible presence of impurities from the manufacturing process.

2. Materials and Methods

2.1. Experimental

2.1.1. Materials

Polycrystalline lovastatin powder (>99%) was purchased from Molekula, Durham, UK and used without further processing. The manufacturer indicates ~1 wt. % impurities, of which 0.5 wt. % are a related compound obtained as a side product from the synthesis.

2.1.2. Infra-Red Spectroscopy

IR spectra were collected using a Thermo Fischer Nicolet FTIR (UK) in attenuated total reflectance (ATR) mode through a diamond crystal (650–4000 cm⁻¹, 4 cm⁻¹ resolution, 32 scans integrated per spectrum).

2.1.3. FT–Raman Spectroscopy

FT–Raman spectra were collected using a Bruker MultiRam (Coventry, UK) equipped with a 1064 nm laser and a liquid nitrogen-cooled Ge detector (30–3600 cm⁻¹, 4 cm⁻¹ resolution, 4 scans integrated per spectrum). The spectra were normalized and 7-point FFT smoothing applied to remove high-frequency thermal noise from the detector.

2.1.4. X-Ray Photoelectron Spectroscopy (XPS)

XPS measurements were performed with a SPECS EnviroESCA (Scanwel, Bala, Gwynedd, UK) using a monochromated Al K_α X-ray source operating at 42 W and a Phoibos 150 NAP analyzer with delay-line detector. The as-received powder was lightly pressed into a pellet and mounted with adhesive carbon tape onto an SEM stub, which was placed on the sample stage. The sample was positioned at a working distance of ~300 µm from the nozzle to the transfer lens stage and was illuminated with a beam spot size of ~300 µm. Gas-phase charge compensation was achieved with 5 mbar Ar in the analysis chamber during measurements. Survey spectra were collected in 1 scan with 100 eV pass energy, 0.5 eV step size and 0.1 s dwell time, while high-resolution core-level spectra were collected in 4 scans with 20 eV pass energy, 0.1 eV step size and 0.4 s dwell time. Peak positions were charge-corrected in binding energy by aligning the aliphatic contribution to the C1s spectrum by −2.75 eV to 285.0 eV and applying the same correction to all spectra. Data analysis was done with CasaXPS [41] Version 2.3.22PR1.0. Core-level spectra were fit using a GL(30) line function for all components and a linear background function.

2.2. Computational

2.2.1. Gas-Phase Conformational Analysis and Simulated Spectra

Gas-phase optimization and conformational analysis of the lovastatin molecule was performed using the Gaussian09 code [42] with the B3LYP functional [43] and a 6-31g+(d,p) basis set [44]. This combination of functional and basis set has been reported to be a good compromise between cost and accuracy for conformational analyses [45].

As described in the text, for the conformational analyses, we investigated the three main bond torsions in the molecule, viz. the “alkyl” C-C torsion around the arm with the pyranyl ring (τ₁), rotation of the O-H bond to different orientations (τ₂), and the “ester” C-O torsion about the ester arm (S-butanoate group) (τ₃). The torsion angles about these three groups were varied over a series of fixed angles in 10° steps, and all other degrees of freedom were optimized at each point. The resulting potential energy surface (PES) from each scan was then inspected to identify the global conformational minimum and three additional low-energy local minima.

The Kohn–Sham orbitals were optimized with the default tolerances of 10⁻⁶ and 10⁻⁸ a.u. on the maximum and root-mean-square (RMS) changes in the density matrix. Geometry optimizations were carried out to the default tolerances of 4.5 × 10⁻⁴ and 3 × 10⁻⁴ a.u. on the maximum and RMS force, and 1.8 × 10⁻³ and 1.2 × 10⁻³ a.u. on the maximum and RMS displacements. The default integration grids were used.

Frequency calculations were performed on the four energy minima to confirm them to be stationary points and to compute theoretical IR and Raman spectra. The Natural Atomic Orbitals (NAO) analysis was used to estimate the core atomic energy levels $E_{\mathrm{NAO}}$, which were taken as binding energies to calculate theoretical XPS spectra [46], i.e.,

$$BE=-E_{\mathrm{NAO}}$$

We note that this is a rather crude approximation—it assumes that $E$ = 0 is the vacuum level and that the photoelectron is completely ejected from the molecule, and does not take into account, e.g., electronic relaxation to stabilize the core hole. In practice, matching the binding energies to measured spectra required offsets of ~11.5 and 18 eV for the C 1s and O 1s binding energies. Simulated spectra were obtained by convolving the individual NAO contribution with a Gaussian function with a full width at half maxima (FWHM) of 0.2 eV to emulate typical instrument broadening, and then summing the resulting series of peak functions.

The electrostatic potentials and electron densities of the conformers were generated using Gaussian09 and visualized using the VESTA software [47]. The Multiwfn software [48] was used to study the reduced density gradients (RDGs), generated on high-quality grids with ~2 $\times$ 10⁶ points, to identify key intramolecular interactions. The isosurfaces were visualized using the Visual Molecular Dynamics (VMD) software [49] in conjunction with scripts provided with Multiwfn.

2.2.2. Solid-State Molecular Dynamics and Simulated Spectra

Periodic density functional theory (DFT) calculations were carried out using the pseudopotential plane wave approach implemented in the Vienna Ab Initio Simulation Package (VASP) code [50].

Calculations were carried out using the PBE generalized-gradient approximation (GGA) functional [51] with the DFT-D3 dispersion correction [52] (i.e., PBE+D3). The electronic structure was modeled using a plane wave basis set with a kinetic energy cutoff of 800 eV and a $\mathsf{\Gamma }$-centered Monkhorst–Pack k-point grid [53] with 3 $\times$ 1 $\times$ 1 subdivisions. These settings were found through explicit convergence testing to converge the total energy and external pressure to < 1 meV atom⁻¹ and 1 kbar (0.1 GPa), respectively. Projector augmented wave (PAW) pseudopotentials [54,55] were used to model the ion cores, with the H 1s and C/O 2s and 2p electrons included in the valence shell.

The starting point for our calculations was the most recent experimental crystal structure with the disorder components removed (CCDC identifier: CEKBEZ01) [2], which was fully optimized to tolerances of 10⁻⁸ eV on the total energy and 10⁻² eV Å⁻¹ on the forces. With these chosen tolerances, optimization with PBE+D3 did not change the cell volume or lattice constants compared to the starting experimental structure.

To study the disorder in the solid, constant-volume Born–Oppenheimer molecular dynamics (BOMD) simulations were performed at 300 K in the single unit cell and a 2 $\times$ 1 $\times$ 1 supercell expansion, for 100 and 25 ps, respectively. The expansion is along the short lattice vector and was chosen as a balance between system size and minimizing artefacts from periodic boundary conditions in the dynamics. The simulation temperature was controlled using a Nosé–Hoover thermostat [56] with the default coupling frequency of 40 timesteps. To speed up the simulations, the plane wave cutoff was reduced to 600 eV, the k-point sampling was reduced to the $\mathsf{\Gamma }$ point, and the mass of the H atoms was increased to 2 (i.e., switching H for D) in order to use a larger MD timestep of 1 fs. Based on explicit convergence tests, the reduced cutoff and k-point sampling produce errors of 0.2 eV molec⁻¹ (~3 meV atom⁻¹) in the total energy and 5.6 kbar in the pressure, but lead to a negligible change in the structure, as determined by comparing the pair distribution functions.

Lattice dynamics calculations were performed on the optimized structure to obtain the $\mathsf{\Gamma }$-point phonon frequencies using the supercell finite displacement method [57] implemented in the Phonopy code [58] and a step size of 10⁻² Å. The absence of imaginary modes in these calculations confirmed that the molecules in the optimized crystal structures are at an energy minimum. The infrared and Raman activities of the modes were then computed from the eigenvectors using the Phonopy-Spectroscopy package [59]. The Born effective charge tensors and dielectric constants required to compute the IR and Raman activities were obtained using the density functional perturbation theory (DFPT) routines in VASP [60]. For the Raman activities, the dielectric constant derivatives were computed using numerical differentiation with a two-point finite difference stencil and a step size of 10⁻²$\times$ the eigenvector norm. Simulated IR and Raman spectra were generated as sums of Lorentzian peak functions using the calculated frequencies and activities and a nominal full width at half maximum (FWHM) of 10 cm⁻¹.

During the geometry optimization and MD simulations, the PAW projection was performed in real space and the precision of the charge density grids was set automatically to avoid aliasing errors. For the MD simulation, the precision was reduced slightly to speed up the calculations. For the lattice dynamics and DFPT calculations, the PAW projection was performed in reciprocal space to ensure accurate forces, and an additional support grid with 8$\times$ as many points was used to reconstruct the augmentation charges.

3. Results

3.1. Conformational Flexibility

3.1.1. Gas-Phase Conformational Analysis

To investigate the conformational flexibility in lovastatin, we first carried out a set of gas-phase calculations to determine the potential energy surfaces (PESs) associated with three major rotatable bonds in the molecule (Figure 1a). A 2D version of the molecular representation can be found in Figure 1e.

As can be seen in Figure 1b–d, the PES scans identify a global minimum and three new low-energy local minima. Two of the new minima were found from the alkyl torsion scan and the third from the hydroxy group orientation scan, while the ester torsion did not show any additional low-energy conformers. The four minima are identified as Conformer 1–4 on the PES scans and labeled in the order of their relative energies from lowest to highest. While we believe that the ab initio conformational analysis has identified a possible global minimum, multi-dimensional PES scan or a gas-phase MD analysis might reveal others. We believe that the present structure, being closest to that of the bulk crystal structure, however, is sufficient for the present study. The relative energies $\Delta E$ and torsion angles τ₁–τ₃ for each of the conformers are given together with other selected geometric parameters in

Table 1. Relative energies, occurrence probabilities at $T$ = 290 K, dihedrals τ₁–τ₃ and selected bond lengths and angles for the four low-energy gas-phase conformers of lovastatin identified in Figure 1.

	Conformer 1	Conformer 2	Conformer 3	Conformer 4
ΔE (kJ mol−1)	0	2.54	4.54	5.70
Occurrence probability pn at T = 290 K	0.627	0.218	0.095	0.059
τ1 (°)	177	180	177	92.1
τ2 (°)	−61.9	60.5	62.2	−61.3
τ3 (°)	0.4	−2.0	−2.9	3.8
$\mathsf{\Phi }$(C19-C28-C20-C17) (butyl) (°)	66.0	66.1	66.0	66.8
$\mathsf{\Phi }$(C10-C28-C20-C17) (butyl, hexahydronaphthenyl) (°)	−172.3	−167.3	−167.6	−165.1
$\mathsf{\Phi }$(O48-C11-C5-C4) (between S-butanoate/pyrnyl groups) (°)	90.7	91.2	90.3	92.1
d(C=O) (S-butanoate) (Å)	1.217	1.217	1.217	1.216
d(C-O) (S-butanoate) (Å)	1.352	1.351	1.352	1.357
d(C=O) (pyranyl) (Å)	1.211	1.210	1.211	1.212
d(C-O) (pyranyl) (Å)	1.353	1.357	1.353	1.351
d(O-H) (pyranyl) (Å)	0.966	0.967	0.966	0.966

. The PES scans also identified three additional high-energy local minima, denoted by asterisks in Figure 1b–d, but these are sufficiently high in energy that we would not expect them to be appreciably populated at room temperature. The PES shown in Figure 1d also indicates a high energy barrier for the rotation of the butanoate group, as previously reported [34].

While we found no significant differences in the bond lengths between the four conformers, each has distinctly different dihedral angles. The alkyl torsion τ₁ represents the angle that the pyranyl group makes with the hexahydronaphathenyl group. The global minimum Conformer 1 has a near-identical dihedral to Conformer 3, and the major difference between them is the angle that the butyl group makes with the hexahydronaphththenyl group ($\mathsf{\Phi }$(C10-C28-C20-C17)). The τ₁ dihedrals of Conformers 2 and 3 differ by only 2° but the hydroxyl group is rotated by 180° (τ₂). In Conformer 4, the τ₁ dihedral is ~90° smaller than the other conformers, which results in the S-butanoate and the pyranyl groups running parallel to each other. An overlay of the structures of the four low-energy conformers and a molecule extracted from the experimental crystal structure can be found in Figure S1 (Supplementary Materials).

Given the relative energies $\Delta E_n$ of the $n$ = 4 conformers, the occurrence probability $P_n$ at a temperature $T$ can be obtained from the expression:

$$P_n\left(T\right)=\frac{1}{Z\left(T\right)}\exp \left[-\Delta E_n/k_{\mathrm{B}}T\right]=\frac{\exp \left[-\Delta E_n/k_{\mathrm{B}}T\right]}{{{\displaystyle \sum }}_n\exp \left[-\Delta E_n/k_{\mathrm{B}}T\right]}$$

where $k_{\mathrm{B}}$ is the Boltzmann constant and $Z\left(T\right)$ is the thermodynamic partition function. At 298 K, this analysis predicts that 62. % of a sample of isolated lovastatin molecules would adopt Conformer 1, 21.8% would adopt Conformer 2, which is 2.54 kJ mol⁻¹ higher in energy, and a much smaller 9.5 and 5.9% would adopt Conformers 3 and 4 due to their higher relative energies of 4.54 and 5.70 kJ mol⁻¹, respectively. We would therefore expect the main flexibility in the molecule to be around the orientation of the OH group and the secondary butyl group. Previous studies [34] using inelastic neutron scattering and solid-state NMR suggested that the main contributions to the spectroscopic features are from methyl dynamics, however, but were unable to resolve the OH dynamics. As shown in the PES scans in Figure 1, the occurrence probabilities of the higher-energy conformers are very strongly temperature-dependent. Excepting the OH group rotation τ₂, the bond torsion listed in Table 1 also does not suggest any angle change greater than 120°, which explains why lovastatin does not show conformational polymorphism [9].

3.1.2. Solid-State Molecular Dynamics

To investigate the conformational flexibility giving rise to the dynamic disorder in bulk lovastatin, we first performed short MD simulations at 300 K on the lovastatin unit cell and a 2 $\times$ 1 $\times$ 1 supercell expansion (Figure 2). Plotting the total energy, temperature, cell pressure and root-mean-square displacement (RMSD) over the two trajectories confirms that the 100 and 25 ps simulation times for the two models were sufficient to reach equilibration (Figures S2–S5).

The static (optimized) structure and the dynamic structures during the MD simulations were compared by computing the pair distribution functions (PDFs) $g\left(r\right)$ (Figure 2a). The $g\left(r\right)$ are calculated for a given structure according to [61]:

$$g\left(r\right)=\frac{1}{4\pi \rho \Delta r}\frac{1}{N}{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^N\frac{\delta \left(r-r_{ij}\right)}{r_{ij}^2}$$

where $N$ is the number of atoms, $\rho =N/V$ is the atomic number density, the double sum runs over atomic pairs $i$ and $j$ separated by a distance $r_{ij}$ and the bond distances are accumulated into a histogram with a bin width $\Delta r$. The $g\left(r\right)$ expresses the probability of finding an atom at a distance $r$ from a reference atom relative to that expected for a homogenous distribution at the average atomic density. Peaks in the $g\left(r\right)$ thus indicate preferred interatomic distances.

In the optimized lovastatin crystal, the complex molecular structure results in sharp features at 1.1, 1.5, 1.8, 2.2 and 2.6 Å, corresponding to five distinct neighbor distances. Comparing the PDFs of the optimized structure to those averaged over the two MD simulations appears to show only thermal broadening, and does not provide any direct evidence for significant conformational changes that might show up, e.g., as shoulder features. Comparison of PDFs averaged over shorter time windows again confirms rapid structural equilibration from the initial model.

We further examined the thermal motion by computing histograms of the atomic positions projected along the crystallographic b and c axes over the two trajectories (Figure 2b,c). Comparison of the two plots suggests a considerably larger range of thermal motion in the trajectory run with the larger supercell, despite there being little evidence for changes in the local structure in the PDFs. The supercell is expanded along the short crystallographic a axis, so we infer from this that the motion of adjacent molecules along this direction is correlated such that movement of a molecule in one unit cell facilitates movement of molecules in the neighboring cells.

We also measured a set of dihedral angles and bond lengths describing the conformation of a representative molecule in each trajectory over the simulation, using the Visual Molecular Dynamics (VMD) software [49]. We measured the τ₁–τ₃ dihedrals investigated in the gas-phase conformational analyses (c.f. Figure 1a) together with an angle quantifying the orientation of the butyl group on the S-butanoate moiety. We also tracked the O-H$\cdots$O H-bond length and an angle defining the O-H group orientation. These data are shown in Figures S6–S8.

While we do observe some thermal variation in the τ₁–τ₃ dihedrals during the simulations, they oscillate about an average value. There is, however, clear evidence for the butyl group adopting two distinct conformations over the two trajectories, suggesting that this is the most flexible part of the molecule in the solid state, as suggested by inelastic neutron scattering data [34]. We observe two conformational changes in the 100 ps trajectory of the single unit cell, and one conformational change during the 25 ps trajectory of the supercell. The higher frequency of changes in the larger supercell may again be taken as evidence that these conformational changes are enabled by cooperative movements of molecules along the a direction, in keeping with the histogram plots in Figure 2b,c. However, conformational changes are inherently stochastic, so this conclusion should be treated with caution given the short simulation times.

Based on the gas-phase calculations, we would expect also to see changes in the orientation of the O-H group. However, our analysis of the MD trajectories shows that the O-H group has a strong preference for remaining in one orientation, with changes in the torsion angle and the O-H$\cdots$O H-bond distance suggest regular, unsuccessful attempts at reaching the rotated conformer. This is likely due to the intermolecular H-bond formed in the solid state, which is not present in the gas-phase models, adjusting the PES for the rotation. We also used the H-bond analysis routine in VMD to count the number of H-bond interactions in the simulation cell using a distance cutoff of 3 Å (Figure S9). This shows that the smaller single-cell model generally has one H-bond per molecule, i.e., all bonds are saturated, whereas the larger supercell model generally has slightly less than one 1 H-bond per molecule. This suggests that, as for the butyl conformation, the dynamics of molecules in neighboring unit cells along the shorter a direction may enable some additional conformational flexibility around the O-H group.

It is, of course, also possible that our simulation is simply not long enough to observe conformational changes. Assuming that a conformational change can be modeled as a first-order reaction, the number of events $n$ per unit time can be estimated from the Arrhenius equation as:

$$n=A\exp \left[-E_A/RT\right]$$

(1)

where $A$ is the attempt frequency, $E_A$ is the activation energy and $R$ is the gas constant. The gas-phase conformational analyses discussed in the following section suggest energy barriers of between 6 kJ mol⁻¹ for rotation of the OH group to 60 kJ mol⁻¹ for the rotation of the S-butanoate about the C-O bond (c.f. Figure 1). If one assumes that the conformational changes are driven by vibrations with a “ballpark” frequency of 1000 cm⁻¹ (30 THz), the corresponding values of $n$ would range from 2.7 to 10⁻⁹ ps⁻¹. We might therefore expect to see the low-energy OH conformational change many times during the simulation, but we would need to simulate for at least nine orders of magnitude longer (i.e., milliseconds) to see the S-butanoate conformation change, if this motion is indeed possible in the solid state. The fact that we do not see evidence for complete OH rotation, in either of the supercells, again strongly suggests that the activation energy may be significantly higher in the solid state than in the gas phase.

3.1.3. Intramolecular Interactions

Lovastatin is a large molecule with a range of possible torsions in the side chains, which can result in different types of intramolecular interactions in different conformers.

Analysis of the electrostatic potentials (ESPs) of the low-energy conformers identified from the gas-phase PES scans shows that all four show subtle electrostatic interactions between the S-butanoate C-O oxygen and the C-OH oxygen on the pyranyl group with H atoms on adjacent C atoms (Figure S10).

We also studied the non-covalent interactions (NCI) in the four conformers from the sign of the reduced density gradient (RDG; Figure 3). The isosurfaces in Figure 3 are colored based on the sign of the second derivative of the electron density (the Laplacian) such that blue indicates H-bonding interactions (-ve Laplacian), green indicates van der Waals (vdW) interactions (Laplacian close to zero), and red indicates steric interactions (+ve Laplacian). The four conformers show an increase in vdW interactions (indicated by green isosurfaces) between the S-butanoate and pyranyl groups with increasing conformer energy. This relates directly to the smaller τ₁ dihedral in Conformer 4, as noted above. There are also increased vdW interactions between the butyl group of the S-butanoate and the hexahydronaphathaenyl group in Conformer 4.

Taken together, these results show that the conformational PES of lovastatin is governed mainly by vdW interactions and does not have a strong electrostatic interaction that preferentially stabilizes the higher-energy conformations or results in any particular intramolecular configurations being “locked in”, as is observed in some organic molecules [62].

3.2. Vibrational Spectroscopy

The P2₁2₁2₁ space group has the 222 (D₂) point group with four irreducible representations (irreps), viz. A, B₁, B₂ and B₃ (Table S1) [63]. Of these, all four may be Raman-active, and the three B irreps may also be IR-active. Ignoring the zero-frequency acoustic modes, the four molecules in the unit cell result in 3$n_a$ − 3 = 777 vibrations, which span an irreducible representation of 195 A + 194 B₁ + 194 B₂ + 194 B₃ at the $\mathsf{\Gamma }$-point. Therefore, as previously reported [34], we anticipate up to 582 IR-active vibrations and 777 Raman-active vibrations, which suggests that both types of spectroscopy could potentially provide rich information about the conformations of molecules in the bulk crystal and at the crystal surfaces, as well as information on functional group terminations at the crystal surfaces.

3.2.1. IR Spectroscopy

The IR spectrum collected from lovastatin powder is shown in Figure 4 together with the solid-state spectrum from lattice dynamics modeling and the simulated spectra of each of the four low-energy conformers identified in the gas-phase PES scans.

In general, our measurements broadly agree with those reported in the literature [36], and a full set of tentative assignments [64] can be found in

Table 2. Tentative IR peak assignments for the reference experimental spectrum and the simulated spectra of solid-state lovastatin and the four low-energy gas-phase conformers identified from the PES scan. Abbreviations: sh—shoulder, w—weak, vw—very weak.

Solid-State	Experimental	Tentative Assignment	Conformer 1 (Gas Phase)	Conformer 2 (Gas Phase)	Conformer 3 (Gas Phase)	Conformer 4 (Gas Phase)
-	655	O-H out-of-plane deformation	652	656	-	-
790 (br, w), 730 (br, w)	698, 714 (w), 739, 752, 760, 782, 796 (sh), 802	C-H out-of-plane deformation in conjugated vinylenes	808, 784, 756, 744, 704	808, 772, 756, 744, 708	808, 796, 772, 756, 744, 708	812, 756, 748, 692
870 (br, w)	833 (w), 841 (w), 870	Aliphatic C-CH3 stretches	876 (sh), 836	-	876	872
976 (w), 966(w)	890, 903, 923, 939 (w), 950 (w), 968	Out-of-plane C-H deformation in cis vinylenes	960, 916, 896 (sh), 884	976 (w), 960 (w), 928, 912, 888,	976, 960, 920 (w), 888	924, 900, 884
1063, 1051, 1031	990, 1014, 1028, 1033 (sh), 1043, 1054	Aliphatic C-H deformation	1156, 1140, 1128, 1112, 1096, 1076, 1068, 1056, 1040, 1020, 1004, 984	1144, 1136, 1124 (w), 1120, 1096, 1084, 1068 (w), 1056, 1036, 1032, 1020 (w), 1004, 984	1116, 1096, 1084, 1068 (sh), 1056 (sh), 1036, 1020, 1004, 988, 984	1108, 1092, 1072, 1056, 1036, 1004, 992
1149, 1115, 1105	1074, 1080, 1107, 1111, 1115, 1126, 1153	Coupled C-O stretches	1340 (s), 1368 (w)	1376 (s), 1340 (br),	1320, 1336	1340 (br), 1372 (w)
1210 (sh), 1199	1168, 1193, 1215	Asymmetric C-O-C stretching	1224 (vw), 1204 (s), 1180 (sh), 1176, 1168 (sh)	1224 (w), 1194 (s), 1176, 1160 (w)	1224 (w), 1200 (w), 1172, 1160 (sh)	1224 (s), 1200 (w), 1180, 1172 (sh)
1285 (br), 1264 (br, sh), 1246 (br)	1246, 1261, 1297, 1329, 1360, 1369	In-plane O-H deformation and symmetric C-H deformation	1272, 1252, 1296,	1320, 1292, 1276, 1260, 1244 (sh)	1276, 1252, 1296,	1248 (sh), 1268, 1256, 1292, 1316,
1395 (w), 1442 (br)	1382, 1444	In-plane O-H deformation coupled with C-H wagging	1420 (sh), 1412, 1384, 1392	1452, 1428, 1424, 1412 (sh), 1384 (br, sh)	1424 (sh), 1412, 1384	1424 (sh), 1412, 1384
	1455, 1461	Asymmetric C-H deformation	1516, 1500, 1484	1520 (sh), 1508, 1500, 1488 (sh)	1508, 1500, 1488 (sh)	1516, 1500, 1488
1664	1698	C=C symmetric stretching in dienes	1677, 1706 (vw)	1677, 1706 (vw)	1677, 1707 (vw)	1675, 1704 (vw)
1697, 1677 (sh)	1711, 1724	C=O in esters and secondary amides (also seen in thiol esters)	1808, 1776	1812, 1772	1808, 1772	1800, 1772
-	2173 (w)	C-H overtones	-	-	-	-
-	2323 (w), 2343 (w), 2362 (w)	Free S-H stretching?	-	-	-	-
2972, 2942, 2918	2989 (sh), 2966, 2950, 2930, 2864	Aliphatic C-H stretching (symmetric and asymmetric in CH2 and CH3)	3052 (sh), 3036, 3012, 2992	3048 (sh), 3036, 3012, 2988	3052 (sh), 3036, 3012, 2988	3052 (s), 3040,3012, 2988
3079, 3041 (s, br)	3014, 3048 (w)	Olefinic C-H stretching	3172, 3132, 3116, 3112 (sh), 3104, 3096 (sh), 3080	3172, 3136, 3120, 3112 (sh), 3100, 3084 (w), 3072 (sh)	3172, 3132, 3120, 3104, 3080	3172, 3144, 3128, 3116, 3104, 3088, 3072 (w)
3538	3538	O-H stretching (with intermolecular H-bonding)	3824	3820	3824	3824

.

The solid-state spectrum from lattice dynamics calculations is a good match with the experimental spectrum, especially at low wavenumbers, although there is a notable shift in the bands in the C-H and C=O regions. These wavenumber-dependent shifts may be due either to a systematic error in the predicted frequencies, or to the fact that the single optimized geometry that we performed the calculations on does not accurately reflect the positional disorder in the crystal.

The calculated vibrational frequencies of the gas-phase conformers are generally higher than those obtained from the solid-state calculations, which may be ascribed to the intermolecular interactions present in the solid-state model but absent in the gas-phase calculations—for example, the observed H-bonding and the CH$\cdots$HC and CH$\cdots$O interactions [34]. In the measured spectrum, the pyranyl group C=O vibration occurs at 1724 cm⁻¹, 10 cm⁻¹ higher than the S-butanoate C=O due to the effect of the ring strain. These vibrations are predicted at 50–70 cm⁻¹ higher wavenumbers in the gas-phase calculations, which we again ascribe to the absence of intermolecular H-bonding. Similarly, coupled C-O stretches at 1500 cm⁻¹ in the gas-phase models appear at much lower wavenumbers between 1050 and 1150 cm⁻¹ in the solid state, again due to constraints imposed by intermolecular H-bonding. The C=C symmetric stretching associated with the diene functionality on the hexahydronaphthenyl group appears in the experimental spectrum as an intense band at 1698 cm⁻¹. The solid-state calculation predicts a similar intensity pattern, while the gas-phase spectra predict a very low intensity, which suggests that interactions between hexahydronaphthenyl groups on adjacent lovastatin molecules in the solid state give rise to appreciable changes in dipole moment. The presence of very weak C-H overtone bands at ~2173 cm⁻¹, together with the presence of vdW interactions suggested by the NCI analysis, is further confirmation for the CH…HC interactions previously suggested using INS studies [34].

Finally, weak features at 2323, 2343, 2362 cm⁻¹ may suggest the presence of impurities with free thiol groups (Figure S12). This is consistent with the ~0.5 wt. % of impurities indicated by the manufacturer, which could include the S-H-containing malonyl-CoA molecule [37]. The acid form of lovastatin is also known to sometimes be present as an impurity [65], but we see no evidence for the presence of carboxylic acid monomers or dimers in our measured spectra. There could also be a contribution to the C=C symmetric stretch region of the spectrum from C=N vibrations of in malonyl-CoA impurities.

Having confirmed good agreement between the measurements and the simulated IR spectrum of bulk lovastatin, we proceeded to compare our measurements to the calculated spectra of the gas-phase conformers, to identify features that may signify differences in the conformation of surface molecules.

The (110) and (020) crystal faces of lovastatin could potentially expose the O-H and olefinic groups of the hexahydronaphthenyl moiety [2], depending on the choice of solvent and preparation method. The solid-state calculations predict a sharp O-H stretch at 3538 cm⁻¹, whereas the corresponding feature in the O-H stretch region of the experimental spectrum is broader and less intense. This can be ascribed to the variable H-bond strength between the O-H with C=O bond on the S-butanoate group in bulk, as well as to surface groups in potentially different chemical environments. Of the four gas-phase conformers, Conformers 1, 3 and 4 have the O-H group in a similar orientation to molecules in bulk, whereas in Conformer 2, the bond is rotated by 180 °. This is predicted to shift the O-H peak by 4 cm⁻¹ in the gas phase compared to the other conformers, despite a negligible 10⁻³ Å change in the bond length. Since the O-H feature in the experimental spectrum is considerably broader than this difference, we cannot discount the presence of Conformer 2 on the surface based on this peak position alone. Furthermore, the intensity patterns of the peaks in the C-H stretching region between 3000 and 3200 cm⁻¹, which contains olefinic C-H and terminal CH₃ stretches, are notably different between the measured and simulated solid-state spectrum, even taking into account the nominal line broadening (10 cm⁻¹) used in the latter. The intensity pattern in the experimental spectrum is, however, very similar to those predicted for the gas-phase Conformer 2 (marked by the green asterisk inside the red box in Figure 4a), which supports appreciable amounts of this conformer being present at the surface.

The lower-wavenumber C-H stretches at 2988 cm⁻¹ and 2992 cm⁻¹ in the gas-phase conformers correspond to C-H stretches on the hexahydronaphathenyl group (C10 in Figure 1e), suggesting that hexahydronaphthenyl groups may be exposed at the surface, as indicated by the olefinic C-H stretches as well. These bands are present in the experimental spectrum but not in the calculated solid-state spectrum, which suggests that hexahydronaphthenyl groups may be exposed at the surface.

The frequency and intensity pattern of the in-plane O-H and symmetric C-H deformations, which occur at 1246, 1261, 1297, 1329, 1360 and 1369 cm⁻¹ in both the experimental and simulated solid-state spectra, are also very similar to the gas-phase spectrum of Conformer 2. The experimental spectra also show O-H out-of-plane deformations at 655 cm⁻¹, which are predicted at 656 cm⁻¹ in the Conformer 2 spectrum but at 652 cm⁻¹ in Conformer 1. This provides further evidence for the possible presence of Conformer 2 on the exposed crystal surfaces, or at least indicates that we cannot discount this possibility based on the spectroscopy.

3.2.2. Raman Spectroscopy

The Raman spectrum of lovastatin powder is shown in Figure 5 alongside the calculated solid-state spectrum and the calculated spectra of the four low-energy gas-phase conformers. A full set of tentative assignments [64] can be found in

Table 3. Tentative Raman peak assignments for the reference experimental spectrum (1064 nm) and the simulated spectra of solid-state lovastatin and the four low-energy gas-phase conformers identified from the PES scan. Abbreviations: br—broad, sh—shoulder, s—sharp, w—weak, vw—very weak.

Solid-State	Experimental	Tentative Assignment	Conformer 1 (Gas Phase)	Conformer 2 (Gas Phase)	Conformer 3 (Gas Phase)	Conformer 4 (Gas Phase)
185, 138, 84, 76	173, 142, 110, 73	Lattice phonon modes	Phonon bands not observed in gas-phase molecules
224, 268, 297		(Unsaturated) CH3 torsional vibration	300, 292, 268	296, 284, 268	300, 280, 264	300, 280, 260
343 (sh), 347, 374, 399	404, 380, 357	C-O out-of-plane deformation vibration	376, 368, 356, 340, 316 (vw)	364, 356, 344,	368, 356, 340,	360, 348, 316
570, 536, 527, 490 (sh), 480, 454, 432	Not visible above FT noise	Aliphatic C-CH3 stretches/ring deformation/C-O in-plane deformation vibration	584, 548, 536, 512, 488, 472, 432, 400	588, 560, 548, 528, 516, 508, 484, 472, 436, 404	588, 560, 548, 528, 508, 484, 472, 436, 404	616, 592, 560, 532, 516, 484, 464, 428, 416, 400
631	655	O-H out-of-plane deformation	652	640	640	656
785, 757, 735, 710, 687, 664	703	C-H out-of-plane deformation in conjugated vinylenes	788, 772, 756, 744, 720, 704, 652	796, 772, 756 (w), 740, 720 (sh), 708, 664	796, 772, 756, 744, 720, 708, 664	796, 756, 748, 720, 692
868, 829, 818	873, 831	Aliphatic C-CH3 stretches	852, 840, 804 (sh), 800	844 (sh), 836 (sh), 832, 808 (sh)	848 (sh), 836, 808	860,856, 840, 812
926, 900, 884	Not visible above FT noise	Out-of-plane C-H deformation in cis vinylenes	988, 984 (sh), 960, 948, 932, 916, 884	984, 976, 960, 948 (w), 932 (w), 928 (w), 912, 896, 888	976 (sh), 960, 932 (w), 912 (w), 888 (w)	972, 956, 940, 920, 900, 884, 872
1076, 1062, 1037, 1029, 1022, 976, 957,	990, 1014, 1028, 1033 (sh), 1043, 1054	Aliphatic C-H deformation	1096, 1076 (sh), 1068, 1056, 1044, 1004	1096, 1080, 1072, 1056, 1044, 1032 (sh), 1004, 992	1096, 1080, 1068, 1056, 1044, 1032, 1004, 992	1092, 1084, 1072 1048, 1032, 992
1151, 1115	Not visible above FT noise	Coupled C-O stretches	1188, 1168, 1156, 1140, 1128, 1112, 1108 (all w)	1188, 1176, 1160, 1144, 1132, 1108 (all w)	1188, 1172, 1160, 1144, 1140, 1132, 1116, 1108 (all w)	1188, 1172, 1160, 1156 (sh), 1144, 1132, 1128, 1108 (all w)
1241, 1217, 1200, 1189	Not visible above FT noise	Asymmetric C-O-C stretching	1260 (sh), 1256, 1248 (sh), 1224, 1220 (w), 1212, 1204	1260 (sh), 1256, 1244, 1228, 1224, 1216, 1204, 1196 (vw)	1260, 1252, 1224, 1204	1256, 1244, 1236. 1216, 1200
1338 (w), 1320 (w), 1298 (sh), 1285	Not visible above FT noise	In-plane O-H deformation and symmetric C-H deformation	1340, 1324, 1304 (w), 1292 (w), 1276 (sh), 1272	1344, 1336, 1328, 1324, 1304 (w), 1276 (sh), 1272	1344, 1340, 1326, 1304 (w), 1276, 1272	1336, 1324, 1316 (w), 1268
1391, 1366 (w, br)	1377, 1404	In-plane O-H deformation coupled with C-H wagging	1384, 1376, 1360	1384, 1376, 1368, 1364	1384, 1364	1372, 1352
1462 (sh), 1436, 1416 (sh)	1446, 1467	Asymmetric C-H deformation	1516 (sh), 1500, 1484 (sh), 1456, 1436, 1432, 1412, 1400, 1396	1516 (br, sh), 1500 (br), 1484 (sh), 1452, 1436, 1432, 1396	1516, 1500 (br), 1496, 1460, 1436, 1432, 1412, 1396	1516 (sh), 1508 (sh), 1500 (br) 1488, 1460, 1432, 1428 (sh), 1416, 1392
1641, 1663 (sh), 1676 (sh)	1647 (br)	C=C symmetric stretching in dienes, C=O in esters?	1808, 1716 (sh), 1712 (sh), 1708, 1704 (s), 1700 (sh)	1812, 1716 (sh), 1712 (sh), 1708 (s), 1704, 1700 (sh)	1808, 1712 (sh), 1708, 1704, 1700 (sh)	1800, 1708 (sh), 1704, 1700 (sh)
3008 (sh), 2981, 2964, 2944 (sh), 2917	2989 (sh), 2966, 2950, 2930, 2864	Aliphatic C-H stretching (both symmetric and asymmetric In CH2 and CH3)	3060 (sh), 3048, 3040, 3012, 2988	3060 (sh), 3052, 3044, 3036, 3012, 2988	3052, 3046, 3012, 2988	3060 (sh), 3052, 3040, 3028 (sh), 3012, 2988
3082, 3063, 3045	3015 (s)	Olefinic C-H stretching	3172, 3148, 3136 (sh), 3132, 3116, 3104, 3080 (w)	3172, 3148 (w), 3136, 3120, 3100, 3084 (w)	3172, 3148, 3136, 3120, 3104, 3076 (w)	3172, 3148, 3132, 3104, 3088 (sh)
3538	Not visible above FT noise	O-H stretching (with intermolecular H-bonding)	3828	3820	3828	3828

.

As observed in the IR spectra, the calculated solid-state spectrum broadly agrees with the experimental measurements, but with the C-H vibrations and C=O stretches in the latter occurring at lower and higher wavenumbers, respectively. The calculated Raman spectra of the gas-phase conformers again highlight that the vibrations in the gas phase generally occur at higher frequencies than in the solid state.

As in the IR spectra, the simulated solid-state Raman spectrum shows a clear O-H peak at 3538 cm⁻¹ which is not seen in the measured spectrum (Figure 5a). This discrepancy could be due to our choice of a 1064 nm laser—since the Raman intensity scales as the fourth power of the wavelength used [66], it is possible that there was insufficient photon flux to excite this particular vibration. Notably, the intensity pattern of the C-H stretches predicted from the solid-state calculations is similar to that of Conformer 1, whereas the intensity pattern in the measured spectrum resembles a mixture of Conformers 1 and 2. The bulk conformation is close to Conformer 1, albeit with a small conformational adjustment [9], whereas our MD simulations suggested that the reoriented O-H bond in Conformer 2 is likely to be strongly disfavored in the bulk. This observation is therefore again consistent with Conformer 2 being present at the surface, as suggested by the IR spectra.

The C=C stretching peaks occur at 1647 cm⁻¹ in the experimental spectrum (Figure 5b). The gas-phase calculations predict that the C=O stretching at ~1808 cm⁻¹ has much lower intensity than the C=C stretching. This, together with the fact that the experimental data were collected using a 1064 nm laser, could suggest that the C=O stretches were perhaps not observed experimentally. In the calculated solid-state spectra, the C=O stretches appear as broad shoulders at 1663 and 1676 cm⁻¹.

The thermal background noise in the spectra prohibits reliable identification of the thiol groups observed in the IR measurements. Features in the fingerprint region from 1100 to 500 cm⁻¹ also had a poor signal-to-noise ratio and could only be reliably discerned after applying 7-point FFT smoothing to de-noise the spectra. The intensity patterns generally resemble the calculated solid-state spectrum (Figure 5c). However, the smoothing of already low-intensity features means that the peak positions may not be reliable, so we refrain from making assignments. The intensity pattern of the C-O out-of-plane deformations between 300 and 400 cm⁻¹ is similar to that predicted in the simulated solid-state spectrum and in the simulated gas-phase spectrum of Conformer 1.

The prominent phonon peak at 110 cm⁻¹ in the experimental spectrum (Figure 5d) is notably absent in the solid-state calculations. There are several possible origins for this discrepancy. The first is that PBE-D3 may incorrectly predict the change in polarizability along the low-frequency modes. The second is that the calculations predict too high an intensity for the strongest bands, so that the relative intensities of the phonon bands are too weak. A third possibility is that the positional disorder in the real structure and consequent symmetry breaking could allow for phonon modes away from the Brillouin zone center to become Raman-active, resulting in an increase in the intensity that is not accounted for in our calculations on a single unit cell. The phonon region of the lovastatin vibrational spectrum has previously been investigated using terahertz spectroscopy [34] and low-wavenumber Raman [28], and these studies found a generally good match between experiment and theory. Our calculations are generally consistent with these previous studies, with calculated vibrational modes at or close to features in the measured terahertz and inelastic neutron scattering spectra in Ref. [34].

Overall, the Raman measurements do not provide a conclusive identification of the conformer(s) present at the surfaces, but do support the inference from the IR measurements that appreciable amounts of Conformer 2 may be present at the surface. Polarized Raman experiments on single crystals with a shorter-wavelength laser might be used in conjunction with modeling, as in our previous study [27], to better probe the functional group orientations at the crystal surfaces. However, since lovastatin tends to form needle-like crystals [67], growing single crystals with sufficient surface area for a typical beam spot size is likely to prove challenging.

3.2.3. X-ray Photoelectron Spectroscopy (XPS)

Given its high surface sensitivity, we also investigated whether XPS measurements might be useful for investigating surface conformations. As described in the Materials and Methods section, the C1s and O1s spectra of the four low-energy conformers identified from the gas-phase PES scans were calculated by treating the atomic energies obtained using the method of Natural Atomic Orbitals as binding energies. We stress once again that this is a very approximate method, and that binding energy shifts of 11.5 and 18 eV needed to be applied to the C1s and O1s spectra, respectively, to match the experimental energy ranges. We therefore use the calculations for qualitative comparison only.

Figure 6a,b show how the simulated C1s and O1s spectra of the lowest-energy gas-phase conformer, Conformer 1, are obtained by summing contributions from the individual atoms. The 11.5 eV shift adjusts the spectra so that the aliphatic C1s binding energy is centered at a reference value of 285.0 eV, similar to the calibration method that we have used in previous XPS studies of organic materials [68]. Similarly, the 18 eV shift applied to the O1s spectra adjusts the C-O peak to match the typical binding energy of 534.5 eV.

After shifting, the C1s component from the C=C in Conformer 1 is predicted to occur at around 284.5 eV, and the C1s peaks arising from aliphatic carbons are predicted to span a range of 284.8–285.4 eV, centered at 285 eV (Figure 6a). The C-C-O C atoms are predicted to give rise to C1s peaks at 285.8 eV, which is a 0.4 eV higher binding energy than the other aliphatic carbons. The C1s peaks arising from the C-O atoms are predicted to produce two distinct peaks at 286.8 eV and 287.4 eV, while peaks arising from O-C=O are centered at 289.0 eV. In both cases where the C is bonded to an O, the peaks occur at a higher binding energy for atoms in the pyranyl ring compared to those in the S-butanoate group. The simulated O 1s spectrum of Conformer 1 (Figure 6b) shows features from the O-C=O atoms centered at 531.5 eV, while those due to C-O atoms are centered at 534.5 eV.

The simulated C1s and O1s spectra (Figure 6c,d) of Conformers 1 and 3 are extremely similar, while the C1s features from the C-C-O atoms in Conformer 2 are predicted to shift to a 0.3 eV lower binding energy. However, the resulting changes to the overall C1s spectrum, indicated by the black boxes in Figure 6c,d, are of the order of 0.1–0.3 eV and may not be detectable. On the other hand, Conformer 4 shows a clear shift in the C1s and O1s features arising from the C-O and O-C=O atoms of the S-butanoate group, which can be related to the increased vdW interactions identified in the analysis in Figure 3, together with a subtle increase in the C-O bond length (c.f. Table 1). These shifts are of the order of 0.5 eV or greater and hence likely could be detected experimentally if an appreciable population of Conformer 4 were to be present.

The survey NAP XP spectrum of lovastatin (Figure 7a) indicates an elemental composition of 74.6% C and 11.5% O with 6.6% N and 7.3% Ar from residual air and charge compensation gas, respectively. Based on the lovastatin stoichiometry, the elemental composition is expected to be 83% C and 17% O. Assuming that the entire O signal in the survey spectrum can be attributed to lovastatin, there would be ~25 at. % excess C present. A common source of such excess C is adventitious carbon contamination from the environment, which is commonly found [68] on materials that have been exposed to air. There may also be impurities in the as-received material, such as N in malonyl CoA, as noted previously.

The C1s core-level spectrum of lovastatin was analyzed by peak fitting using four components, without constraints, reflecting the expected contributions of the four primary C moieties in the molecule to the total peak area. These are as follows: 4 sp² C atoms in the in the two C=C bonds (17%), 15 sp³ aliphatic carbon atoms (62%), 3 C (12%) in a single bond to a heteroatom (alcohol and ester) and 2 C atoms in the C=O carbonyl groups (8%). The components have similar FWHM (~1.8 eV) except for the component reflecting C atoms in the ester and alcohol functional groups, which has a slightly narrower FWHM (1.3 eV). The narrower FWHM suggests that these C atoms are in more similar chemical environments than the C atoms contributing to other components of the spectrum. The fit approximates the relative theoretical contributions of the four clusters of functional groups plus those from adventitious carbon to the spectral envelope. The fit also broadly agrees with the trends in the calculated C 1s spectra of the gas-phase lovastatin conformers (c.f. Figure 6c).

The O1s core-level spectrum was fitted with 2 peaks, with an equal FWHM broader than those of the C 1s components, reflecting 3 O atoms at lower BE contributing 60% to the total intensity and 2 O atoms at higher BE contributing to the remaining 40% of the overall intensity. Based on the molecular structure of lovastatin, one expects the contributions to be reversed, with the 2 O atoms from the carbonyl group contributing to the spectral envelope at a lower binding energy than the 3 O atoms in the ester and alcohol species. This result would also be in line with the relative contributions of these species observed in the calculated O1s spectra. However, neither of these interpretations account for the impact of non-covalent inter- or intramolecular interactions on the electron density of O atoms of lovastatin in the solid state. For example, computational studies of crystalline lovastatin [2] have identified that non-covalent intermolecular interactions involving the O atoms in the carbonyl and alcohol functional groups make significant contributions to stabilizing the crystal lattice. These include a short intermolecular contact between the carbonyl and an acidic proton on adjacent pyranyl rings and a long H-bond between a carbonyl group on the butanoate ester and an alcohol on an adjacent pyranyl ring. Redistribution of electron density at the O atoms from these electrostatic interactions is well known to shift the binding energies of the core-level electrons [69] by a similar magnitude to the shifts observed in the O1s spectrum presented here. The relatively broader peaks in the O1s spectrum also suggest that the 5 O species might be distributed over more than just two defined states. Resolving these details is, however, impossible without a better systematic understanding of core-level binding energy shifts as a function of the exact proton positions within the hydrogen bonds. The influence of these non-covalent interactions on the core-level binding energies is stronger than the effects of conformational change (Figure 6), and a deeper understanding of the dynamic and static disorder of protons, as well as the depths of the PES of H-bonding, will be needed to identify the preferred conformation of lovastatin in the solid state.

4. Discussion

Gas-phase conformational analyses of the main bond torsions identify four low-energy conformers that would be expected to show appreciable populations at room temperature. Analysis of the electrostatic potentials and non-covalent interactions does not suggest any particular electrostatic interactions that would make one conformer lower in energy than the other, so the energy differences likely arise from a subtle balance of intramolecular interactions. The relative energies suggest a 21.8% occupation of a second conformer at 290 K, with the O-H group on the pyranyl ring pointing 180° to that in the other conformers. Across the four low-energy gas-phase conformers, the two main chains of the molecule comprising the S-butanoate and the lactone do not exhibit changes to any torsion that are greater than 120°, thereby concurring with previous studies [9] about the absence of conformational polymorphism. Solid-state molecular dynamics simulations, on the other hand, suggest that the O-H group is locked in position in bulk by a strong H-bonding interaction with the neighboring C=O group. The disorder observed in the crystal structure is instead due to the butyl group on the S-butanoate moiety, and the MD simulations suggest co-operative movement molecules in adjacent unit cells along the short a direction.

A comparison of measured vibrational spectra to simulated spectra of bulk lovastatin and the gas-phase conformations suggests that an appreciable amount of the second conformer may be found at crystal surfaces, where the group is able to rotate freely. In particular, the olefinic C-H stretching vibrations from the hexahydronaphathenyl group, the coupled in-plane O-H deformation and symmetric C-H deformation and the appearance of the out-of-place O-H deformation at 656 cm⁻¹ in the experimental IR spectrum support the presence of the second conformer at the surface. The intensity pattern of the C-H region in the FT–Raman spectrum also suggests a combination of the two lowest-energy gas-phase conformers, again supporting the presence of the second conformer at the surface. The implication of this is that lovastatin crystal surfaces likely expose free O-H groups and are therefore not chemically inert.

While experimental and computed XPS generally agree, we found little observable difference between the two prominent low-energy gas-phase conformers. Better distinction between conformers and impurities present on the true surface and sub-surface will require more detailed studies using variable photon energy and angle-dependent XPS experiments [70].

We also observed evidence in our measurements for the presence of malonyl CoA impurities in our commercial lovastatin sample. XPS measurements indicate a small at. % of N, which is difficult to identify in IR spectra due to the relatively weak N-H vibrations, while S-H vibrations are evident in the IR spectra but are not observed in XPS due to the lower sensitivity to S. These two measurements are therefore complementary.

Overall, the combination of gas-phase and solid-state modeling and spectroscopy has provided valuable insight into the conformational flexibility of the lovastatin molecule in bulk and at the particle surfaces, in a powder sample. In the future, our protocol could be improved on by performing polarized Raman measurements, as well as angle-resolved XPS measurements on clean single-crystal surfaces. This combined approach could be applied to a wide range of molecular solids, and we therefore hope that our work will serve as a first step towards a robust protocol for the routine analysis of particle surfaces.