(623f) Describing Polarization in Hybrid Ab Initio/Empirical Force-Field Models for Crystal Structure Prediction

Prediction of polymorphism is an important problem in solid-state chemistry due to the dependence of crystal physico-chemical properties on the underlying crystal structure. In service of this, Crystal Structure Prediction (CSP) seeks to predict the putative polymorphs of a solid system and rank them by their relative thermodynamic stability. CSP methods which combine ab initio calculations with empirical force-fields have proved to be effective at determining crystal thermodynamic stability [1, 2]. However, such approaches natively neglect the effects of polarization, which are induced by the electric field generated from molecules in the crystalline lattice. Modelling polarization is crucial for the accurate description of intermolecular interactions like hydrogen bonds [3, 4], and its omission significantly hampers model accuracy. In this work, we investigate approaches to incorporate polarization into the hybrid ab initio/empirical force-field (HAIEFF) CSP model, at different levels of fidelity.

HAIEFF models rely on isolated-molecule ab initio computations conducted in the gas-phase to describe, among other things, electrostatic interactions in the crystal. The isolated nature of these calculations makes them cheap to evaluate, enabling the efficient use of HAIEFF models in repeated energy minimizations and in studies involving multicomponent crystals. However, a consequence of using isolated-molecule calculations is the neglect of polarization, which can limit the accuracy of HAIEFF models. To embed polarization within this framework, it is desirable to instead perform these isolated-molecule computations with approximations of the crystalline electric field. To that end, two induction models are studied in this work. An implicit representation of the crystalline electric field can be attained using a Polarizable Continuum Model (PCM) [5], and this has been employed within CSP in the past [1, 6]. To achieve a more-detailed representation of the crystalline electric field, an electronic-embedding method termed the self-consistent electronic response to point-charges (SCERP) model [3, 4] has also been proposed. To our knowledge, this latter model has not been actively applied in CSP. For meaningful appraisal of these induction models, compatible empirical force-fields have been independently parametrized as part of this work. Failing to carry out such parametrization would entail incongruent use of the HAIEFF model, which introduces indeterminant errors into the modelling results. By fitting an additional force-field in the absence of either polarization method, we also consider a third model (denoted FF_ind), wherein the effects of polarization are completely borne by the force-field parameters. These three models (FF_ind, PCM, and SCERP) are benchmarked through two CSP case studies.

The efficacy of the three models is first analyzed through their performance in the force-field parametrization. The parameter estimation training data comprises geometry and energy data determined from high-accuracy periodic density functional theory (DFT) calculations [7]. This data is available for a large and diverse set of 604 organic molecular crystals, with and without intermolecular hydrogen bonds. Formulating the HAIEFF force-field as a Buckingham potential was found to give suboptimal performance, particularly in hydrogen bonding crystals. Attempts to accommodate polarization in the force-field led to poorly balanced intermolecular forces and a significant occurrence of molecular fusion. These issues were resolved by instead formulating the force-field as a Mie potential, which is better suited for describing highly attractive interactions like hydrogen bonds. Between the three HAIEFF models, FF_ind had a higher propensity to underbind hydrogen bonding crystals relative to the two models that embed polarization into their ab initio component. Overall, the training data was more accurately reproduced when using PCM or SCERP than when relying solely on the force-field to absorb polarization effects. Between PCM and SCERP, differences in parameter estimation accuracy were marginal.

These three models were then applied in CSP studies of a rigid monohydrate system and a flexible small organic molecule. In both cases, FF_ind performed notably worse than the PCM and SCERP models, with at least one polymorphic form predicted significantly higher in relative energy than is reasonable. Between PCM and SCERP, CSP predictions were generally of similar quality. However, the PCM model failed to capture some subtle differences in polarization between the different crystal geometries, leading to slightly erroneous predictions. In contrast, when using the SCERP model, these nuances were more precisely described, such that the CSP results conformed better to predictions from periodic DFT.

Ultimately, our results reiterate that polarization plays a significant role in determining relative polymorphic stabilities. At the lowest fidelity, the HAIEFF force-field can absorb errors from the absence of a polarization model, however the efficacy of this FF_ind approach is limited. Models like PCM and SCERP that accommodate polarization through the ab initio component exhibit markedly better model efficacy and should be more widely adopted in CSP. In its current state, SCERP provides some marginal improvement over PCM. However, the more physically robust foundations of SCERP offer greater potential for systematic improvement, compared to PCM. Pursuing these avenues may be necessary for HAIEFF models to achieve significant improvement in model accuracy.

References

[1] Reilly, A. M. et al. Report on the sixth blind test of organic crystal structure prediction methods. Acta Crystallographica Section B 2016, 72, 439–459.

[2] Hunnisett, L. M. et al. The seventh blind test of crystal structure prediction: structure generation methods. Acta Crystallographica Section B 2024, 80, 517–547.

[3] Welch, G.; Karamertzanis, P.; Misquitta, A.; Stone, A.; Price, S. Is the Induction Energy Important for Modeling Organic Crystals? Journal of Chemical Theory and Computation – J CHEM THEORY COMPUT 2008, 4, 522–532.

[4] Zhang, Y. Crystal Structure Prediction for Cocrystals and Hydrates: Energy Models and Structure Generation. Ph.D. thesis, Imperial College London, 2023.

[5] Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chemical Reviews 2005, 105, 2999–3094.

[6] Cooper, T. G.; Hejczyk, K. E.; Jones, W.; Day, G. M. Molecular Polarization Effects on the Relative Energies of the Real and Putative Crystal Structures of Valine. Journal of Chemical Theory and Computation 2008, 4, 1795–1805.

[7] Bowskill, D. H.; Tan, B. I.; Keates, A.; Sugden, I. J.; Adjiman, C. S.; Pantelides, C. C. Large-Scale Parameter Estimation for Crystal Structure Prediction. Part 1: Dataset, Methodology, and Implementation. Journal of Chemical Theory and Computation 2024, 20, 10288–10315.