Hybrid classical/machine-learning force fields for the accurate description of molecular condensed-phase systems

Abstract

Electronic structure methods offer in principle accurate predictions of molecular properties, however, their applicability is limited by computational costs. Empirical methods are cheaper, but come with inherent approximations and are dependent on the quality and quantity of training data. The rise of machine learning (ML) force fields (FFs) exacerbates limitations related to training data even further, especially for condensed-phase systems for which the generation of large and high-quality training datasets is difficult. Here, we propose a hybrid ML/classical FF model that is parametrized exclusively on high-quality ab initio data of dimers and monomers in vacuum but is transferable to condensed-phase systems. The proposed hybrid model combines our previous ML-parametrized classical model with ML corrections for situations where classical approximations break down, thus combining the robustness and efficiency of classical FFs with the flexibility of ML. Extensive validation on benchmarking datasets and experimental condensed-phase data, including organic liquids and small-molecule crystal structures, showcases how the proposed approach may promote FF development and unlock the full potential of classical FFs.

Fig. 1. Classification of classical intermolecular interactions. (left) Classical fixed-charge FF with point charges and a Lennard-Jones potential. (middle) Classical polarizable FF such as AMOEBA with the additional inclusion of polarization (∞, 1) and atomic multipoles (2, (0, 1, 2)). (right) Model proposed in this work, which includes a three-body dispersion term (3, 1) and a pairwise ML potential (2, (0, 1, 2)) compared to the polarizable FF. The pairwise ML potential can account for directional interactions in a systematic manner.

Fig. 2. Overview of the ANA2B model. Dotted lines refer to features that depend on the geometry while bold lines to features based on molecular graphs. Blue components refer to intermolecular interactions, red to intramolecular interactions, and grey to shared interactions.

Fig. 3. Results for the lattice energies of systems in the X23 (top panels) and ICE13 datasets (bottom panels) with the ANA2B models. The left column (ANA2B⁰) refers to the model without any treatment of polarization, the middle (ANA2B¹) column shows results for the model that includes only the direct polarization term, and the right column (ANA2B^∞) displays results for the model which includes a full treatment of polarization. Equality ±4.184 kJ mol⁻¹ is indicated by the black lines.

Fig. 4. Results for condensed-phase properties of 57 molecules used in the parametrization and validation of GROMOS 2016H66: density (left) and heat of vaporization (right) for the ANA2B¹ model. Black lines indicate equality ±50 kg m⁻³ and ±4.184 kJ mol⁻¹.

Fig. 5. Stability ranking for the crystal structure for the compounds of the CSP blind tests 3 (VIII, X, XI) and 5 (XVI, XVII, XVIII) using the lattice energy predicted with the ANA2B^∞ model. Each horizontal bar represents the stability of a structure with respect to the most stable structure. Red bars indicate experimental structures. The candidate structures were taken from the corresponding publications.

Fig. 6. Stability ranking for the crystal structure of compound XXII. Each horizontal bar represents the stability of a structure with respect to the most stable structure. The stability is given in kJ per mol per molecule. Candidate structures and rankings for dispersion corrected PBE and PBE0 are taken from ref. . Experimental structures are marked in red.

Fig. 7. Stability ranking for the crystal structure of compound XXIII. Each horizontal bar represents the stability of a structure with respect to the most stable structure. The stability is given in kJ per mol per molecule. Candidate structures and rankings for dispersion corrected PBE and PBE0 are taken from ref. . Experimental structures are marked in color.

Fig. 8. Stability ranking for the crystal structure of compound XXVI. Each horizontal bar represents the stability of a structure with respect to the most stable structure. The stability is given in kJ per mol per molecule. Candidate structures and rankings for dispersion corrected PBE and PBE0 are taken from ref. . Experimental structures are marked in red. Note that PBE0 + MBD is not available for all polymorphs of this structure.