Structure solution using the experimentally determined unit cell constants was performed using a combined molecular dynamics (MD) and density functional theory (DFT) approach. MD simulations were performed using LAMMPS in combination with the CHARMM General Force Field v. 2b7.24,25 Several thousand trial structures were generated by randomly placing one molecule into an expanded unit cell (120% the experimental unit cell volume). During the MD simulation run (70 ps with 1 fs time step), the system was allowed to relax while the unit cell shrunk back to the experimentally determined size. After equilibration at the experimental system size and a final energy minimization step, the atomic positions were extracted. Resultant structures were grouped based on their geometric similarity using the Hausdorff metric. In each structural grouping, characteristic features such as distinct hydrogen bond lengths or π–π interaction distances could be identified. Subsequently, to overcome the shortcomings of classical MD simulations, one representative structure from each group was used as an input for DFT geometry optimization as implemented in the VASP program (version 5.4.1).26–29 The PBE functional was used for exchange and correlation,30 with PAW potentials for all the elements.31 The effects of van der Waals interactions (vdW) were included using the pairwise method of Tkatchenko and Scheffler (TS).32 A plane-wave cutoff energy of 800 eV was employed. Calculations were performed on the primitive cell, containing 30 atoms, and a 3 × 2 × 1 Monkhorst–Pack k-grid to sample the Brillouin zone.33 The total energy during the self-consistency loop of each DFT step was converged to 10–8 eV. Calculations were performed using the experimental volume, relaxing atomic positions down to a threshold of 10–2 eV Å–1 on forces. The final optimized geometry was verified by calculating the phonon dispersion using density functional perturbation theory as implemented in VASP.