Conformation versus cohesion in the relative stabilities of gabapentin polymorphs

Abstract

Polymorphism in molecular solids is driven by the competition between internal molecular conformational strain and external binding forces. The pharmaceutical gabapentin exhibits three polymorphic forms each differing in molecular conformation and intermolecular interactions. To improve the understanding of the polymorphic formation, terahertz spectroscopy and solid-state density functional theory have been applied in the analysis of the factors leading to these particular crystalline arrangements. It was found that the most stable solid-state polymorph contained gabapentin molecules under considerable internal strain. The origins of this unexpected stability are described in terms of intramolecular conformational changes and the amount of cohesive binding energy present in each polymorph. Finally, the temperature dependence of the polymorph stabilities is investigated by calculation of the free–energy curves of each form and determining the solid–solid transition points connecting the various phases.

---

1. Introduction

The formation of distinct arrangements of chemical species in the solid state involves the competition between intermolecular binding forces driving condensation and molecular conformational strain destabilizing the process. The ability to crystallize into multiple solid-state forms, polymorphism,¹ is caused by variations in the balance between these external and internal forces. Therefore, it is possible to have polymorphs with higher relative internal conformational energies becoming the more stable solid-state forms through external cohesive forces countering such energy penalties.^2,3 The tendency to crystallize when burdened with higher conformational strain is exhibited by many compounds, especially larger and more flexible molecules such as pharmaceuticals. A consequence of the interplay of these forces is the creation of metastable forms in which the molecular arrangement in the solid is temporarily stable under certain conditions. These metastable forms are typically kinetically populated, while the stable forms are thermodynamically favored. With time, they will eventually relax into stable forms as the balance of forces tilts in one direction. The conditions to achieve metastability can be related to factors such as pressure or temperature, with the metastable form destabilizing when the factor is modified or removed, since the outside force is what provides the additional energy to resist thermodynamic relaxation. Understanding these energetic details in a polymorphic system can yield insight into all types of polymorphism, including traditional and conformational polymorphism.^1,4–6

Gabapentin^7–10 (Fig. 1), the active ingredient in Neurontin, is prescribed for prevention of seizures in patients with neurological disorders and to relieve neuropathic pain. The polymorphism exhibited by gabapentin provides an excellent example of the balance achieved between the intermolecular and intramolecular forces found in pharmaceutical solids, and the role of metastability in the overall ranking of crystalline stabilities. Gabapentin is known to exist in three zwitterionic anhydrous forms, Form II (P2₁/c, Z = 4, Z′ = 1), Form III (P2₁/c, Z = 4, Z′ = 1), and Form IV (C2/c, Z = 8, Z′ = 1), as shown in Fig. 2, and one hydrate (Form I). Forms II, III, and IV are also known as alpha, beta, and gamma, respectively.⁷ While crystallizing in the same space group, Forms II and III do exhibit conformational polymorphism (Fig. 3), specifically a ring inversion that flips the axial and equatorial positions of the substituents. This significant change in ring orientation provides an opportunity to study polymorph formation as a function of molecular conformation versus intermolecular forces.

Fig. 1 Molecular structure of gabapentin with atomic labeling scheme.

Fig. 2 Unit cell molecular packing of the three polymorphs of gabapentin.

Fig. 3 Side-view comparison of the ring conformations in gabapentin Forms II and III.

Changes in solid-state packing between polymorphs result in variations in the low-frequency vibrations exhibited by the molecules, including both internal (molecular torsions) and external (rotations and translations) motions. Investigation of these low-frequency vibrations can be accomplished using terahertz (far-infrared) spectroscopy, which probes the intermolecular and intramolecular vibrations of molecular solids in the range of approximately 10 to 100 cm⁻¹. The ability to access these vibrations enables terahertz spectroscopy to be used to analyze and understand the polymorphic nature of pharmaceuticals^11–14 or the subtle effects of disorder in a crystalline solid.^15,16 In this study, it has been used to investigate the sub-90 cm⁻¹ vibrations of two polymorphs of gabapentin, Forms II and III, followed by a complete structural and vibrational analysis using quantum mechanical models with periodic boundary conditions.

Solid-state density functional theory (DFT) provides an excellent means for enhancing our understanding of the chemical origins of the observed terahertz spectra by simultaneous simulation of both the intramolecular and intermolecular vibrations that are important in this spectral region. Using solid-state DFT, high correlation with experimental data can be achieved and it is a proven way to interpret and assign terahertz spectra.^17–22 By achieving the complete computational reproduction of an observed terahertz spectrum, as well as the crystal structure of the solid being considered, the solid-state simulations demonstrate that they are able to well represent all of the forces in the sample. These models can then be utilized to obtain the relative energies of the polymorphs in order to rank their overall stabilities, but also more detailed molecular conformational energies and cohesive binding forces can be extracted. Beyond the electronic energies that are manifested as conformation and cohesion, the simulations also allow for the calculation of the thermodynamic properties of each polymorph, such as zero-point energy and entropy, which yield important information about the real-world stabilities of each form.

2. Methods

2.1 Experimental

Gabapentin was purchased from Sigma Aldrich (≥98% purity) and the supplied solid-state form was verified by single-crystal X-ray diffraction to correspond to Form II. Form III was obtained by recrystallization of gabapentin from ethanol at 333 K.^7,10 The preparation of Form IV proved to be difficult and while methods for its recrystallization have been published,^10,23,24 no amount of Form IV could be produced following these procedures. The single-crystal X-ray diffraction data used can be found in the ESI.† For terahertz sample preparation, each polymorph was mixed with a powdered polytetrafluoroethylene (PTFE, Teflon) matrix (Spurlock Specialty Tools, ∼3 micron particles) at a concentration that was dependent on the polymorph. Form II samples used 6.5%, and Form III 1.2% by mass to ensure that measurements were made within an optimal absorption range. The mixtures were then pulverized using a stainless-steel grinder/mixer (Dentsply Rinn 3110-3A) to minimize particle size and reduce both Mie scattering and crystal anisotropy.²⁵ Pressure induced polymorphism is not uncommon in organic molecular crystals²⁶ and polymorphic transformations have been reported for gabapentin¹⁰ after long (15 minutes) grinding periods. No evidence of any transformations was found in the present study after brief (30 seconds) grinding was used. Approximately 0.55 g (total) of the sample mixtures were pressed, at 2000 psi, into pellets with diameters of 13 mm and thicknesses of 2.0 mm. Pure PTFE was pressed into a pellet of equal mass to be used as a blank reference.

The experimental spectra were obtained using a time-domain pulsed terahertz spectrometer based on an amplified Ti:Sapphire femtosecond laser system. Zinc telluride crystals were used for both generation of terahertz radiation by optical rectification and detection by free-space electro-optic sampling.^27,28 A detailed description of the terahertz spectrometer has been reported elsewhere.²⁹ The samples and blank for measurement were held under vacuum in a cryostat with data acquired at both 293 K and 78 K. Samples and blanks were scanned 32 times for each individual data set over a time window of 32 ps (chosen to exclude all pulse reflections) consisting of 3200 data points, which was then symmetrically zero-padded to a total of 6000 data points. The ratio of the Fourier-transformed data sets of the sample and blank resulted in a terahertz spectrum over the range of 10 to 90 cm⁻¹ with a spectral resolution of approximately 1.0 cm⁻¹. Each data set was replicated four times at both temperatures and then averaged to obtain the final spectra reported here. Spectral intensities are reported in units of ε (M⁻¹ cm⁻¹) where molarity is expressed in terms of the concentration of crystallographic unit cells (Z) rather than individual molecules.

2.2 Theoretical

Both solid-state and gas-phase (isolated molecule) calculations were performed for this study. The vibrational analyses, structural simulations, and energetic investigations of the crystalline polymorphs were all completed using solid-state DFT with periodic boundary conditions. While the solid-state calculations were the principal computational method, gas-phase simulations were used as a complementary method to further the understanding of the possible mechanisms underlying the relative polymorphic energies and transformations in the solids.

The solid-state simulations in this work were performed using the CRYSTAL09³⁰ software package utilizing the PBE³¹ density functional in combination with the atom-centered Gaussian-type 6-311G(2d,2p)³² basis set. The total energy convergence criteria were ΔE < 10⁻⁸ hartree for geometry optimizations and ΔE < 10⁻¹¹ hartree for frequency calculations. All structural optimizations were performed without limits on atomic positions or unit cell dimensions, other than those imposed by space group symmetry, and were begun using starting structures obtained by experimental X-ray diffraction measurements. The Form II and III structures were measured for this work at 90 K using single crystal X-ray diffraction. Because of the difficulty in crystallization of Form IV, the previously reported X-ray structure (at 173 K) was used as the starting point for its simulation.⁷ The radial and angular distributions for DFT integration were defined by a pruned (75, 974) grid. Truncation tolerances for Coulomb and HF exchange integrals were defined as 10⁻⁷, 10⁻⁷, 10⁻⁷, 10⁻⁷, 10⁻¹⁴ hartree (TOLINTEG command^30,33). A shrinking factor of 6 (64 k points in the irreducible Brillouin zone) was determined after sampling and monitoring of the total energy convergence as a function of k-point count in reciprocal space according to the Pack–Monkhorst method.³⁴ Normal-mode frequencies and infrared intensities were then calculated for the optimized structures. The frequency of each normal mode was calculated within the harmonic approximation by numerical differentiation of the analytical gradient of the potential energy with respect to atomic position.³³ The infrared intensities for each normal mode were calculated from the dipole moment derivatives (dμ/dQ) determined using the Berry phase technique of calculating Born charges as polarization differences between equilibrium and distorted geometries.^30,35 Mode descriptions and assignments were made by visual inspection of the eigenvector atomic displacements for each normal mode.

The explicit consideration of weak non-covalent interactions in molecular solids can be a significant factor in achieving valid simulations of crystalline structure and dynamics.^36,37 London-type dispersion force corrections are important augmentations to include in solid-state DFT investigations since these forces are generally underestimated in commonly used density functionals, but play an integral role in solid-state characteristics. The solid-state DFT approach used in this study has been supplemented with corrections for London-type dispersion forces using a semi-empirical method proposed by Grimme,³⁸ and then later modified for use in the solid-state by Civalleri et al.³⁶ A global scaling factor (s₆) of 0.60 was used for all three polymorphs, and was obtained through comparison of the calculated unit cell parameters (a, b, c, α, β, γ, and volume) and the experimental 90 K X-ray data taken for Forms II and III. Only comparisons to low-temperature structures were made due to the sensitivity of the scaling factor to experimental temperature.³⁷

Gas-phase calculations were completed with the Gaussian09³⁹ software package, utilizing the M06-2X⁴⁰ density functional and the same basis set as used in the solid-state simulations. This particular functional was chosen based upon its performance versus other functionals as reported by Zhao and Truhlar.⁴⁰ Additionally, in order to check the quality of some of the DFT results, gas-phase Møller–Plesset perturbation theory (MP2)⁴¹ was also applied for treatment of electron correlation. All gas-phase calculations used pre-optimized molecular structures extracted from the solid-state simulations as starting points. All geometry optimizations allowed the molecules to fully relax with no symmetry restrictions, while the single-point energy calculations held the molecules in rigid conformations. The default total energy convergence criteria for both DFT and MP2 simulations was used (ΔE < 10⁻⁶ hartree) for geometry optimizations and single-point energy calculations. The DFT grid was set to program option “ultrafine” for all simulations using M06-2X. A polarizable continuum solvent model⁴² was also utilized to further ease the comparison of the calculations to experiment. Intermediate transition states were also considered using the STQN method^43,44 to locate and simulate the transition state structure between two different molecular conformations of gabapentin, with the QST2 command from Gaussian09. The starting and ending structures of the transition state search were selected as corresponding to the Form II and Form III conformations, and the program was allowed to find the intermediate transition state.

3. Results and discussion

3.1 Terahertz spectroscopy

Terahertz spectroscopy directly probes the low-frequency intermolecular and intramolecular vibrations of molecules in the solid-state offering insight into the low-energy interactions between molecules. The analytical power of terahertz spectroscopy is illustrated here through the observed differences in the spectra of gabapentin polymorphs, enabling unambiguous identification of the different forms. The 78 K and 293 K terahertz spectra from 10 to 90 cm⁻¹ of Forms II and III of gabapentin are shown in Fig. 4 and the observed vibrational frequencies are listed in Table 1. In addition to the clear spectral features seen in the spectrum of Form II, there is likely a strong absorption just outside the spectrometer range, resulting in the sharply rising baseline beyond ∼80 cm⁻¹. The terahertz spectrum of gabapentin Form IV was not obtained because of the aforementioned difficulty in making this polymorph. The low-temperature spectra exhibit sharpened features due to the decrease in the populated vibrational states of the molecules in the solid-state, making the 78 K spectra essential for the proper assignment of the vibrational modes. All of the 78 K spectral features for gabapentin Form II exhibit a slight shift to higher energy versus their room-temperature positions, a common event in the terahertz spectroscopy of cooled solids.^11,29,45 In contrast, one of the 78 K spectral features observed in Form III (33.33 cm⁻¹) appears at a lower energy than found at room temperature. This unusual phenomenon has been observed before in the terahertz spectroscopy of sucrose⁴⁶ and irbesartan,¹⁶ but it is rare. In order to better understand the nature of the spectral features and to study the origins of the energetic differences, solid-state DFT has been employed in the spectral analyses.

Fig. 4 Terahertz spectra of gabapentin crystallized in two different forms. 293 K data is shown in gray and 78 K data is shown in black.

Table 1 Observed terahertz vibrational frequencies (cm⁻¹) of gabapentin Forms II and III at 293 K and 78 K and the approximate correlations between the two temperatures

Form II		Form III
293 K	78 K	293 K	78 K
39.44	42.22	35.55	33.33
48.88	52.22	41.66	43.88
81.66	86.10	56.66	58.88
		64.99	64.99
		70.55	77.77
		86.66	87.77

---

3.2 Solid-State DFT simulations

3.2.1 Structural analyses. X-ray data at 90 K was only available for two of the three gabapentin polymorphs (Forms II and III); therefore the quality of the structural reproduction is focused on these. The same parameters used for the simulations of Forms II and III were used in the simulation of Form IV, with the structure being compared to the 173 K X-ray data found in the literature.⁷ The difference in X-ray diffraction measurement temperatures is unfortunate, as the temperature variation will add some small level of error into the structural comparison, but unavoidable due to the difficulty in growing Form IV. Structure simulations were based on full optimizations (using PBE/6-311 G(2d,2p)), without limits on the unit cell dimensions. Root mean squared deviations (RMSDs) of bond lengths, bond angles, and dihedral angles (Table 2) were used to determine the quality of the reproduced structures. The RMSD values represented here exclude all hydrogen atoms. The quality of the unit cell parameters (a, b, c, and volume) was also evaluated and the results are listed in Table 3. The structural reproduction of Form II is the best of the three polymorphs, with the lowest RMSD values and unit cell parameters only differing by 0.45% from experiment. Form III also optimized well, with similar RMSD values when compared to Form II for bond lengths and bond angles. However, the torsion angle errors for Form III are more than double that of Form II. The optimized unit cell parameters for Form III also deviated from experiment by more than a factor of two as compared to the Form II results, but the 1.12% error still indicates that it is an overall good structural reproduction. While the calculated structures of Forms II and III were of fine quality, the simulation of Form IV was not as successful. The RMSD values for Form IV are the worst of the three polymorphs, and in combination with the 6.12% change in the unit cell parameters, indicate that Form IV did not optimize to a structure in good correlation with the reported X-ray data.⁷ This result may be interpreted as a failure of the applied theory, but it is unexpected that the theory would work well for only two of the three polymorphs given their similarities. The computational results may be indicating that there are some issues with the reported structure or that there may be a temperature dependence in the Form IV stability (the simulations are effectively 0 K). These possibilities can be evaluated by investigating the energy components that contribute to the relative stabilities of the polymorphs and also the vibrational motions of these forms as the means for determining their thermodynamic parameters (entropy, free energy, etc.). The relative polymorphic stability information deepens our understanding of the polymorphic formation mechanisms, while also providing possible explanations for the difficulty in obtaining a good structural reproduction of Form IV.

Table 2 Evaluation of the structural reproduction using PBE/6-311 G(2d,2p) by root mean squared deviations (RMSDs) of each polymorph of gabapentin, excluding hydrogens. Bond lengths are in angstroms (Å), while bond and dihedral angles are in degrees (°)

	Bonds	Angles	Dihedrals
Form II	0.0050	0.318	0.679
Form III	0.0080	0.365	1.555
Form IV	0.0151	0.637	1.769

---

Table 3 Comparison of the 90 K experimental unit cell dimensions and the simulated unit cell dimensions of Forms II and III. Form IV is compared to 173 K experimental data.⁷ All simulations were completed using the global scaling factor (s₆) of 0.60

	Form II		Form III		Form IV
	Experimental	Theoretical	Experimental	Theoretical	Experimental	Theoretical
a (Å)	5.86	5.92	14.51	14.86	30.55	33.84
b (Å)	6.92	6.92	6.63	6.63	5.93	6.00
c (Å)	22.25	22.18	9.81	9.75	10.88	11.00
α (deg)	90.00	90.00	90.00	90.00	90.00	90.00
β (deg)	90.00	89.48	105.73	106.40	108.32	111.07
γ (deg)	90.00	90.00	90.00	90.00	90.00	90.00
V (Å^3)	903.50	907.79	908.87	921.32	1870.58	2083.90
% error in Unit cell parameters		0.45		1.12		6.12

---

3.2.2 Vibrational analyses. The simulated terahertz spectra of Forms II and III are displayed in Fig. 5. The calculated vibrational modes (frequencies and intensities for modes ≤90 cm⁻¹), mode descriptions, and mode assignments can be found in Table 4. The intensities were calculated in units of km mol⁻¹, but were converted into spectral absorption units (M⁻¹ cm⁻¹) for direct comparison (no arbitrary scalars applied) with experimental data. This was achieved by using an empirically-determined Lorentzian line shape with a full-width half-maximum (fwhm) value of 2.37 cm⁻¹ for Form II and 3.44 cm⁻¹ for Form III.

Fig. 5 Overlay of the 78 K experimental terahertz spectra of gabapentin Forms II and III (dot) and the spectra simulated by PBE/6-311 G(2d,2p) (solid).

Table 4 Vibrational frequencies (cm⁻¹) and intensities (km mol⁻¹) of the calculated IR-active vibrational modes for gabapentin Forms II (α) and III (β), below 100 cm^−1a

Form II					Form III
Mode	Freq	Int	Mode description	Obs	Mode	Freq	Int	Mode description	Obs
a Freq: theoretical frequency, int: mode intensity, obs: observed experimental frequency, Rj: external rotation about j-axis, Tj: external translational along j-axis, IM: internal motion, *Tentatively assigned to rising absorption at end of spectral bandwidth.
1α	26.93	0.01	Rb	—	1β	32.16	2.90	Rc	33.33
2α	40.33	0.03	Tc	—	2β	40.52	6.38	Rb	43.88
3α	41.14	0.14	Tb	42.22	3β	56.48	1.07	Rb	—
4α	52.25	1.59	Rb	52.22	4β	59.27	8.15	Tb	58.88
5α	87.64	0.20	Ra	86.10	5β	60.76	14.95	Rc	64.99
6α	90.81	0.01	Rc	—	6β	74.52	40.20	50% Rb, 50% IM	77.77
7α	102.49	6.73	IM	*	7β	85.12	4.50	Ra	87.77
					8β	87.42	14.96	Ra

---

The experimental terahertz absorptions in the Form II spectrum were assigned using the calculated frequencies of Table 4, but only the vibrations exhibiting significant IR intensities were considered. The lowest-frequency experimental feature at 42.22 cm⁻¹ corresponds to 3α, the vibrational mode calculated at 41.14 cm⁻¹, which originates from a translation along the crystallographic b-axis. The calculated vibrational mode 4α, at 52.25 cm⁻¹, is produced by an external rotational motion around the b-axis, and is assigned to the experimental mode found at 52.22 cm⁻¹. The highest-frequency experimental spectral feature at 86.10 cm⁻¹ originates from the vibrational mode 5α calculated at 87.64 cm⁻¹, which is produced by an external rotation about the a-axis.

All of the Form III experimental spectral features could be readily correlated with the calculated vibrational modes listed in Table 4. The lowest frequency mode at 33.33 cm⁻¹ can be attributed to the vibrational mode 1β calculated at 32.16 cm⁻¹, generated by an external rotation around the c-axis. The calculated vibrational mode 2β at 40.52 cm⁻¹ was generated by an external rotation about the crystallographic b-axis, and is assigned to the experimental spectral feature at 43.88 cm⁻¹. The experimental mode at 58.88 cm⁻¹ can be attributed to the calculated mode 4β, at 59.27 cm⁻¹, which is generated by a rigid-body translation along the b-axis. The calculated mode 5β, at 60.76 cm⁻¹, can be assigned to the spectral feature at 64.99 cm⁻¹, and is produced by an external rotation around the c-axis. The experimental mode at 77.77 cm⁻¹ originates from the vibrational mode 6β calculated at 74.52 cm⁻¹, which is produced by an external rotation about the b-axis coupled with an internal torsion about the C6–C7 (Fig. 1) bond. The highest-frequency experimental spectral feature at 87.77 cm⁻¹ can be assigned to the calculated modes 7β and 8β (at 85.12 and 87.42 cm⁻¹, respectively), both generated by external rotations around the crystallographic a-axis.

The excellent correlation between theory and experiment for both gabapentin polymorphs indicates that the intermolecular forces are well represented by the models. The computational reproduction of both the terahertz spectra and the crystal structures enables the detailed analysis of the polymorph energies to proceed with confidence.

3.2.3 Relative energies of the polymorphs in the solid-state. Probing the origins of the total energies of these molecular crystals by solid-state DFT reveals a coexistence of cohesive and conformational energies that balance each other to achieve the observed forms. A result of this internal conflict is that having lower conformational strain does not always translate to a lower total energy in the solid-state since the influence of cohesive interactions in the solid-state may dominate and vice versa. Examining the differences in the relative energies between the gabapentin polymorphs leads to greater understanding of polymorph formation and stability, especially the ability to overcome energetic penalties caused by molecular conformation.

While solid-state simulations can produce accurate stability rankings,³ finite basis sets impart an error to the final energies, necessitating a correction. In order to generate the best representation of the energy rankings, the basis set superposition error (BSSE) must be estimated and removed from the stability values. The counterpoise method⁴⁷ was used here for approximating the amount of BSSE in the solid-state simulations. While this method can be straightforward for gas-phase systems (limited number of atoms), when dealing with periodic boundary condition simulations, a spatial cutoff must be defined for inclusion of neighboring basis functions in the counterpoise calculation. In the calculation of BSSE, a single molecule was extracted from the already optimized solid-state unit cell and studied using the same theoretical method (PBE/6-311G(2d,2p)) as for the periodic calculations. It was found, through monitoring of the energy differences, that 200 atoms within 5.0 Å of the molecule being evaluated were sufficient limits, and such cut-offs were used in all calculations to gauge BSSE effects. The BSSE (per molecule) was found to be −46.95 kJ mol⁻¹ for Form II, −38.76 kJ mol⁻¹ for Form III, and −42.46 kJ mol⁻¹ for Form IV. While this approach does not recover all of the BSSE, it was found that at least 95% of the BSSE energy is revealed.

The total electronic relative energies (per molecule, Table 5) indicate that Form II was the most stable (0.00 kJ mol⁻¹) polymorph, Form III was next at +1.51 kJ mol⁻¹, and Form IV was the least stable polymorph (+6.51 kJ mol⁻¹). The small energy gap between Forms II and III shows that the two polymorphs are separated by less than ambient room-temperature energy (∼2.4 kJ mol⁻¹ at 293 K). These differences in energy are consistent with Form II being the more common and favored solid-state form of anhydrous gabapentin.⁷ Given the similarities in polymorph energies, contributions from zero-point energy (ZPE) cannot be ignored. Application of the ZPE corrections (obtained from the frequency analyses) does not change the stability ranking, but does change the absolute energy values. The gap between Form II, which is still the most stable polymorph (0.00 kJ mol⁻¹), and Form III decreases by ∼47%, making the two polymorphs very close in energy (now separated by only 0.80 kJ mol⁻¹). Form IV is still the least stable (with an energy of +3.32 kJ mol⁻¹), however, it also is much more stable after inclusion of the ZPE correction.

Table 5 Calculated energies per molecule (kJ mol⁻¹) of three polymorphs of gabapentin, all corrected for BSSE

	Form II	Form III	Form IV
a Zero point energy correction.
Relative total energies	0.00	1.51	6.51
Relative total energies with ZPE^a	0.00	0.80	3.32
Relative conformational energy	90.13	0.00	74.42
Total binding energy	−357.85	−266.22	−335.64

---

The total electronic energy of each polymorph is a result of the conflict between conformational strain and intermolecular binding.^2,16 Therefore, in order to better understand the origins of the relative polymorphic stabilities, both conformational molecular energy and cohesive binding energy have been determined for all three polymorphs. The single-molecule energies were calculated using the same extracted single molecule that was used for BSSE determination, and was studied using the same theoretical method that has been used throughout. Basis set superposition error was also corrected in the cohesive binding energy calculations. Examining the conformational energies reveals that gabapentin molecules in the Form III polymorph have the most stable conformation (0.00 kJ mol⁻¹), with Forms IV and II being considerably less conformationally stable by +74.42 kJ mol⁻¹ and +90.13 kJ mol⁻¹, respectively. Such a large difference in conformational energy between the polymorphs is unexpected, especially after the small differences in the total electronic energy determined in the solid-state. Therefore, while knowing the conformational energies does give some insight into the total energy rankings of the polymorphs, it clearly does not explain everything, resulting in the need for examining the remaining energetic component, cohesive binding.

Intermolecular interactions play a central role in the formation and energetics of polymorphic systems. The large differences between the total energy rankings and the conformational energy rankings for gabapentin indicate that the cohesive binding energies are especially important for understanding polymorph formation in this compound. The binding energy can be determined from the total solid-state energies corrected for BSSE and the energies of the individual isolated molecules extracted from the crystalline solids. The difference between the total solid energy and the summation of its molecular component energies is the cohesive energy of the solid. It was found that Form II has the greatest amount of binding energy per molecule with −357.85 kJ mol⁻¹, while Form III has the least amount of binding energy (−266.22 kJ mol⁻¹). The binding energy of Form IV falls between the other two polymorphs at −335.64 kJ mol⁻¹. These energies represent the collective forces between the molecules in the solid state and include factors such as dipole–dipole interactions, hydrogen bonding, and London dispersion forces. While the total binding energies are the focus of the present work, it is important to note that computational techniques such as Hirshfeld surfaces analysis⁴⁸ could potentially be utilized to further characterize the intermolecular binding energies.

The total solid-state energy of the gabapentin polymorphs derives from a combination of the binding and conformational energies. While the total solid-state energies are in agreement with the experimental polymorphic formation from ethanol⁷ (Form II recrystallizes at room-temperature, Form III at 333 K), the conformational energy seems to indicate that Form III should be the most stable polymorph at room temperature. However, the favorable intermolecular binding energy of Form II, makes it the more stable solid-state polymorph overall. The structural reproductions, solid-state polymorphic energetics, and the experimental difficulty in making Form IV (using the method from Reece and Levindis⁷ produced Form II exclusively), suggest that there is a complex mechanism underlying the formation of gabapentin polymorphs.

3.2.4 Examination of the conformational energy differences. While the relative solid-state energies were valuable for understanding the polymorphic formation in gabapentin, some important aspects of the intramolecular energetics, that are relevant to the overall crystallization process, can be studied using isolated-molecule (gas-phase) simulations, potentially elucidating the polymorph ranking and formation mechanisms. While PBE was used in the solid-state simulations because of its proven reliability, better performing density functionals exist for use in the gas phase, and the M06-2X density functional was chosen for this task.⁴⁰ In order to test the variation in the molecular energies by using the two different density functionals, single-point energy calculations were completed to find the relative conformational energies of the polymorphs, using the geometries of the pre-optimized solid-state molecular conformations. The M06-2X/6-311G(2d,2p) simulations produced results very similar to PBE, with Form III being the most stable conformationally (0.00 kJ mol⁻¹), while Forms IV and II were less stable by +79.85 kJ mol⁻¹ and +92.82 kJ mol⁻¹, respectively. To check the possible influence of explicit electron correlation in the conformational energies, Møller–Plesset perturbation theory (MP2) was used in combination with the 6-311G(2d,2p) basis set to evaluate these same single-point energies. The MP2 results agreed with the conformational energy rankings of both density functionals (Form III: 0.00 kJ mol⁻¹, Form IV: +82.57 kJ mol⁻¹, and Form II: +95.54 kJ mol⁻¹) indicating the reliability of the methods being applied.

While the gas-phase simulations utilizing different theoretical approaches verified the conformational energy rankings, the large energetic gap present between Form III and the other two polymorphs still required additional calculations to be performed in order to understand its origin. The next step beyond isolated-molecule simulations was to include solvent effects that may play a role in the crystallization process. Since all three polymorphs are soluble in ethanol, and all three forms recrystallize from ethanol at different temperatures, the gas-phase calculations were extended using the polarizable continuum solvent model (PCM).⁴² The single-point energy calculations were completed again, with M06-2X/6-311G(2d,2p) and PCM represented ethanol. The relative energy ranking of the forms did not change, however, the absolute values did change significantly, with Form III being the most stable conformationally (0.00 kJ mol⁻¹), Form IV at +38.57 kJ mol⁻¹, and Form II still being the least stable at +47.00 kJ mol⁻¹. The continuum model solvent effectively halved the energy difference between the different gabapentin conformations found in the three polymorphs, suggesting that solvent plays a key role in the polymorph formation mechanisms.

Since zwitterions are generally unstable in the gas-phase,^49–51 the rigid single-point energy simulations were completed first. However, full relaxation of the individual molecules could also provide valuable information to the relative stabilities. The optimizations were performed with M06-2X/6-311G(2d,2p) and PCM ethanol. All of the geometry optimizations lead to a proton transfer from the NH₃⁺ to the COO⁻, turning the zwitterions into their neutral counterparts. While the neutral forms make comparison to the experiment more difficult, the results do provide useful qualitative information. Form III was found to be the most stable at 0.00 kJ mol⁻¹, while Forms II and IV both relaxed to an identical higher-energy structure with an energy +1.58 kJ mol⁻¹ above Form III. These results indicated that only two distinct conformations of gabapentin exist in solution, that of Form II/IV and Form III (Fig. 3). The geometry optimization results also demonstrate that the neighboring molecules in the solid-state are able to hold the gabapentin molecules in unfavorable conformations that easily relax when unrestrained.

Forms II and IV do have different conformations (as is illustrated by the rigid conformational energies), but these structural differences are based solely on the positioning of the COO⁻ and NH₃⁺ groups enforced by crystal packing. When optimized without the constraints of the solid-state structure (intermolecular binding) the two polymorphs arrived at the same lower energy conformation. However, Form III achieved the lowest energy conformation based on its unique ring orientation rather than a small change in substituent interactions. The key difference between the two optimized forms (Form III vs. Form II/IV) is the conformation of the ring, and the only pathway from one polymorph to the other is by a ring inversion. In order to better understand this conformational change, a single-molecule transition state search, using the STQN method, was completed (simulated in ethanol) between the Form III and Form II/IV geometries. The search resulted in a transition state structure consisting of a twisted boat conformation of the ring and an intermediate positioning of the ring substituents. Using this transition structure, the barrier height between the two conformations was determined to be 43.52 kJ mol⁻¹, which is very close to the experimental value for the ring inversion barrier in cyclohexane (∼43.1 kJ mol⁻¹).^52–54 This similarity suggests that the source of the barrier is in the ring itself rather than the attached functional groups.

The use of gas-phase simulations for zwitterionic molecules is less than ideal, because of the difficulty in stabilizing the zwitterion.^49–51 In the case of gabapentin, it was thought that the addition of a discrete ethanol molecule to the calculation, in combination with the PCM solvent, would help prevent the formation of the neutral species. An ethanol molecule was placed between the COO⁻ and NH₃⁺ groups (with the ethanol oxygen near the NH₃⁺ group and the hydroxyl hydrogen near the COO⁻ group) with the intention being to stabilize the zwitterion by blocking the internal proton transfer coordinate. The M06-2X/6-311 G(2d,2p) optimizations allowed the molecules to fully relax and resulted in stable zwitterionic gabapentin molecules. Form III was again found to be the most stable polymorph (0.00 kJ mol⁻¹), while Form II and IV relaxed to the same higher-energy structure (+6.49 kJ mol⁻¹). The same type of STQN transition state search was performed on the structures using the explicit ethanol solvent molecule in combination with the PCM solvent. The barrier height in this case was found to be 43.20 kJ mol⁻¹, essentially unchanged from the neutral gabapentin results, and still in agreement with the experimental value for the ring inversion for cyclohexane. Thus, the zwitterionic or neutral character of gabapentin appears to have little impact on the ring inversion coordinate.

The difference in conformational energy between Forms II and III was found to be ∼90 kJ mol⁻¹ when extracted as rigid molecules from the crystals, with Form III being more stable. However, the total solid-state energy ranking indicates that Form II is more stable, if only by a small margin. The gas-phase calculations reveal the energy barrier found between the two molecular conformations further complicating the formation process of the polymorphs. The large difference in binding energy between the two polymorphs must be sufficient to force a ring inversion, allowing Form II to become the most stable form at room-temperature. Therefore, Form III must only crystallize at higher temperatures, in order to overcome the binding energy differences, as well as the inversion barrier between the conformations. The temperature dependence of the polymorph formation can be seen through an investigation into the thermodynamic properties of the polymorphs, yielding insight into their experimental stabilities.

3.2.5 Thermodynamic properties involved in the polymorphic formation. Results from the solid-state and gas-phase simulations have agreed with the literature⁹ in that Form III is metastable at room temperature. By examining the thermodynamic properties (mainly Gibbs free energy) of the solid-state polymorphs over a series of temperatures, information about this metastability can be gained. The frequency analyses performed in the CRYSTAL09 software package allow for the scanning of the thermodynamic properties versus temperature (ignoring any temperature-induced volume changes) by combining the electronic energy values (zero-point energy corrected) with the thermodynamic quantities (TS, PV, and thermal contribution to the vibrational energy). These temperature scans for each polymorph are shown in Fig. 6 and reveal that at higher temperatures, Form III is the most stable form, while Form II is the least stable. This graph also indicates that there are two solid–solid phase transition points, one from Form II to Form III at ∼85 K and one between Forms II and IV at ∼260 K. Not shown in Fig. 6 is that the Forms III and IV curves eventually cross at a much higher temperature of ∼1850 K. While the absolute values are obviously incorrect when compared to the experimental observations, the temperature ordering of the transitions are in agreement with the previously reported solid–solid polymorphic transformations.⁵⁵ In that study, it was determined that Form II was the most stable at room temperature, Form III formed through a solid–solid transformation from Form II at 358 K, and Form IV was produced by heating Form III to 393 K. Therefore, these analyses can only be viewed qualitatively, and indicate Form II starts as the most stable solid-state polymorph, and as the temperature rises, becomes the least stable polymorph, while Form IV becomes the most stable form at much higher temperatures. In this case, low-temperature differential scanning calorimetry (DSC) could be a useful technique for locating possible transitions below the 20–180 °C range investigated by Lin et al.⁵⁵ The computational thermodynamic study verifies the concept that Form III is metastable,⁹ while also strongly suggesting that Form IV is a metastable form of gabapentin, although it requires much more energy for this metastable state to be significantly populated. The inability to produce the Form IV sample limits the amount of information that can be ascertained about this polymorph, however, the results presented here are supported by the formation of this polymorph from the melt at high temperatures.⁵⁵ The high energy required to make Form IV also explains the difficulty in making the polymorph at lower temperatures in ethanol, since the available energy (via the TS term in the Gibbs free energy) would be insufficient to initiate the Form II/Form IV phase transformation.

Fig. 6 Thermodynamic energy curves of gabapentin Forms II (solid), III (dash), and IV (dot).

4. Conclusions

The terahertz spectra of two gabapentin polymorphs were measured in the spectral range of 10 to 90 cm⁻¹ at room temperature and 78 K. Both polymorphs exhibited distinct spectral features that could be used to uniquely identify each crystalline form. Solid-state DFT simulations allowed for a complete analysis of each terahertz spectrum and enabled assignments to specific solid-state motions. These assignments indicated that the terahertz absorptions are largely due to intermolecular motions of translational and rotational character and therefore highly sensitive to packing variations between the polymorphs.

A combination of gas-phase and solid-state DFT simulations of three anhydrous polymorphs of gabapentin were used to rank their relative stabilities, including partitioning into cohesive binding energy and molecular conformational strain contributions. It was found that while Form II was the most stable solid-state form, it exists in a molecular conformation that is very unstable compared to Form III. This indicates that the intermolecular binding forces contribute considerably to the stability of the polymorphs. Scanning of the Gibbs free energy of each polymorph versus temperature suggests that Form III is the high-temperature stable form, once the differences in binding energy are overcome. Form IV is the most stable form at extremely high temperatures (∼1850 K), but this is beyond the melting point of gabapentin. This investigation into the relative stabilities of the polymorphs of gabapentin suggests that at high temperatures, the conformational energy within the molecules determines the solid-state arrangement, but at low temperatures, intermolecular forces dominate and supercede conformational strain. These results increase our understanding of the complex formation processes of various polymorphs and will help characterize existing and predict new polymorphs in other pharmaceutical solids.

Acknowledgements

This research was funded by a grant from the National Science Foundation CAREER program (CHE-0847405). The authors also thank Syracuse University for continued support.

References

1. T. L. Threlfall, Analyst, 1995, 120, 2435–2460.
2. A. Nangia, Acc. Chem. Res., 2008, 41, 595–604.
3. S. P. Delaney, D. Pan, S. X. Yin, T. M. Smith and T. M. Korter, Cryst. Growth Des., 2013, 13, 2943–2952.
4. E. Weber, Top. Curr. Chem., 1998, 198, 220.
5. J. Bernstein, Polymorphism in Molecular Crystals, Oxford University Press, New York City, U.S., 2002.
6. Polymorphism: In the Pharmaceutical Industry, ed. R. Hilfiker, Wiley-VCH Verlag GmbH & Co. KGaA, 2006.
7. H. A. Reece and D. C. Levendis, Acta Crystallogr., Sect. C: Cryst. Struct. Commun., 2008, 64, o105–o108.
8. J. A. Ibers, Acta Crystallogr., Sect. C: Cryst. Struct. Commun., 2001, 57, 641–643.
9. K. E. Dempah, D. H. Barich, A. M. Kaushal, Z. Zong, S. D. Desai, R. Suryanarayanan, L. Kirsch and E. J. Munson, AAPS PharmSciTech, 2013, 14, 19–28.
10. D. Braga, F. Grepioni, L. Maini, K. Rubini, M. Polito, R. Brescello, L. Cotarca, M. T. Duarte, V. Andre and M. F. M. Piedade, New J. Chem., 2008, 32, 1788–1795.
11. M. D. King, W. D. Buchanan and T. M. Korter, Anal. Chem., 2011, 83, 3786–3792.
12. P. F. Taday, I. V. Bradley, D. D. Arnone and M. Pepper, J. Pharm. Sci., 2003, 92, 831–838.
13. K. Ajito, Y. Ueno, H.-J. Song, E. Tamechika and N. Kukutsu, Mol. Cryst. Liq. Cryst., 2011, 538, 33–38.
14. J. A. Zeitler, P. F. Taday, D. A. Newnham, M. Pepper, K. C. Gordon and T. Rades, J. Pharm. Pharmacol., 2007, 59, 209–223.
15. R. Li, J. A. Zeitler, D. Tomerini, E. P. J. Parrott, L. F. Gladden and G. M. Day, Phys. Chem. Chem. Phys., 2010, 12, 5329–5340.
16. S. P. Delaney, D. Pan, M. Galella, S. X. Yin and T. M. Korter, Cryst. Growth Des., 2012, 12, 5017–5024.
17. P. U. Jepsen and S. J. Clark, Chem. Phys. Lett., 2007, 442, 275–280.
18. M. R. C. Williams, D. J. Aschaffenburg, B. K. Ofori-Okai and C. A. Schmuttenmaer, J. Phys. Chem. B, 2013, 117, 10444–10461.
19. S. P. Delaney, E. M. Witko, T. M. Smith and T. M. Korter, J. Phys. Chem. A, 2012, 116, 8051–8057.
20. E. M. Witko, W. D. Buchanan and T. M. Korter, J. Phys. Chem. A, 2011, 115, 12410–12418.
21. M. D. King, P. M. Hakey and T. M. Korter, J. Phys. Chem. A, 2010, 114, 2945–2953.
22. O. Kambara, K. Tominaga, J.-i. Nishizawa, T. Sasaki, H.-W. Wang and M. Hayashi, Chem. Phys. Lett., 2010, 498, 86–89.
23. Z. Zong, S. D. Desai, A. M. Kaushal, D. H. Barich, H.-S. Huang, E. J. Munson, R. Suryanarayanan and L. E. Kirsch, AAPS PharmSciTech, 2011, 12, 924–931.
24. C.-H. Hsu, W.-T. Ke and S.-Y. Lin, J. Pharm. Pharm. Sci., 2010, 13, 67–77.
25. T. M. Lowry and F. C. Hemmings, J. Soc. Chem. Ind., London, 1920, 39, 101–110.
26. E. Boldyreva, Acta Crystallogr., Sect. A: Found. Crystallogr., 2008, 64, 218–231.
27. A. Rice, Y. Jin, X. F. Ma, X. C. Zhang, D. Bliss, J. Larkin and M. Alexander, Appl. Phys. Lett., 1994, 64, 1324–1326.
28. Q. Wu, M. Litz and X. C. Zhang, Appl. Phys. Lett., 1996, 68, 2924–2926.
29. P. M. Hakey, D. G. Allis, W. Ouellette and T. M. Korter, J. Phys. Chem. A, 2009, 113, 5119–5127.
30. R. Dovesi, V. R. Saunders, C. Roetti, R. Orlando, C. M. Zicovich-Wilson, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. J. Bush, P. D'Arco and M. Llunell, Crystal09 User's Manual, 2009.
31. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
32. R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980, 72, 650–654.
33. F. Pascale, C. M. Zicovich-Wilson, F. L. Gejo, B. Civalleri, R. Orlando and R. Dovesi, J. Comput. Chem., 2004, 25, 888–897.
34. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188–5192.
35. S. Dall'Olio, R. Dovesi and R. Resta, Phys. Rev. B: Condens. Matter, 1997, 56, 10105–10114.
36. B. Civalleri, C. M. Zicovich-Wilson, L. Valenzano and P. Ugliengo, CrystEngComm, 2008, 10, 405–410.
37. M. D. King and T. M. Korter, J. Phys. Chem. A, 2012, 116, 6927–6934.
38. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
39. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, GAUSSIAN 09, Gaussian, Inc., Wallingford, CT, 2009.
40. Y. Zhao and D. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
41. M. Head-Gordon, J. A. Pople and M. J. Frisch, Chem. Phys. Lett., 1988, 153, 503–506.
42. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3093.
43. C. Peng, P. Ayala, H. B. Schlegel and M. J. Frisch, J. Comput. Chem., 1996, 17, 49–56.
44. C. Peng and H. B. Schlegel, Isr. J. Chem., 1994, 33, 449–454.
45. M. D. King, W. D. Buchanan and T. M. Korter, Phys. Chem. Chem. Phys., 2011, 13, 4250–4259.
46. M. Walther, B. M. Fischer and J. P. Uhd, Chem. Phys., 2003, 288, 261–268.
47. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553–566.
48. M. A. Spackman and D. Jayatilaka, CrystEngComm, 2009, 11, 19–32.
49. W. D. Price, R. A. Jockusch and E. R. Williams, J. Am. Chem. Soc., 1997, 119, 11988–11989.
50. R. R. Julian and M. F. Jarrold, J. Phys. Chem. A, 2004, 108, 10861–10864.
51. Y. Ding and K. Krogh-Jespersen, Chem. Phys. Lett., 1992, 199, 261–266.
52. R. K. Harris and N. Sheppard, J. Mol. Spectrosc., 1967, 23, 231–235.
53. R. K. Harris and N. Sheppard, Proc. Chem. Soc., London, 1961, 418–419.
54. F. A. L. Anet, M. Ahmad and L. D. Hall, Proc. Chem. Soc., London, 1964, 145–146.
55. S.-Y. Lin, C.-H. Hsu and W.-T. Ke, Int. J. Pharm., 2010, 396, 83–90.

---

Footnote

† Electronic supplementary information (ESI) available: Single-crystal X-ray diffraction data for Forms II and III taken at 90 K. See DOI: 10.1039/c3ra43887b

---