Capturing Dichotomic Solvent Behavior in Solute-Solvent Reactions with Neural Network Potentials

Abstract

Simulations of chemical reactivity in condensed phase systems represent an ongoing challenge in computational chemistry, where traditional quantum chemical approaches typically struggle with both the size of the system and the potential complexity of the reaction. Here, we introduce a workflow aimed at efficiently training neural network potentials (NNPs) to explore energy barriers in solution at the hybrid density functional theory level. The computational burden associated with training at the PBE0-D3(BJ) level is bypassed through the use of active and transfer learning techniques, whereas extensive sampling of the transition state region is accelerated by well-tempered metadynamics simulations using multiple time step integration. These NNPs serve to explore a puzzling solute-solvent reactivity route involving the ring opening of N-enoxyphthalimide experimentally observed in methanol but not in 2,2,2-trifluoroethanol (TFE). This reaction represents a challenging example characterized by intricate hydrogen bonding networks and structurally ambiguous solvent-sensitive transition states. The methodology successfully delivers detailed free energy surfaces and relative energy barriers in agreement with experiment. These barriers are associated with an ensemble of transition states involving the direct participation of up to five solvent molecules. While this picture contrasts with the single transition state structure assumed by current static models, no drastic qualitative difference is observed between the formed hydrogen bonding networks and the number of participating solvent molecules in methanol or TFE. The dichotomy between the two solvents thus essentially arises from an electronic effect (i.e., distinct nucleophilicity) and from the larger conformational entropy contributions in methanol. This example underscores the critical role that dynamic simulations at the ab initio levels play in capturing the full complexity of solute-solvent interactions.

Scheme 1. (a) Rhodium Catalyzed Reaction of an N-enoxyphthalimide (1) with an Alkene (2) to Form Either Carboamination (3, Red Pathway) or Cyclopropanation Products (4, Blue Pathway). (b) Pathway of Solvent-Mediated Ring Opening of N-enoxyphthalimide via a Stepwise Mechanism

Figure 1

Process for TS structure sampling and database reinforcement. First, direct and baselined GLNNPs are trained on the initial training database (from Section 2.2). Well-tempered metadynamics simulations combined with a multiple-time step integration are then run to explore the TS region. New structures arising from these simulations are categorized as either lying in the extrapolative or interpolative regime by committee analysis. Any structure deemed to be an extrapolation is then extracted, PBE-D3(BJ) energies and forces computed, and the structure is added into the updated training database. This process is repeated iteratively until no structures are deemed to be extrapolations. Once this occurs, a new set of direct GLNNPs (only) is trained and nanosecond-length well-tempered metadynamics simulations run. The Same committee analysis is then used to again check for structures lying in the extrapolative regime. When all structures from these longer simulations are considered to be interpolations, the PBE-D3(BJ) training database is complete and the NNPs can be used to map the free energy surface of the reaction.

Figure 2

(a) Pipeline depicting the transfer learning procedure. Ten percent of structures from the GLNNP database are extracted using the FPS algorithm on the two CVs, s and z, which serve as the initial HLNNP database. NNPs are then trained with initial weights taken from the GLNNPs. (b) Parity plot comparing xTB and PBE0-D3(BJ) energies (in eV) for the initial PBE0-D3(BJ) database in (a). (c) Parity plot comparing PBE-D3(BJ) and PBE0-D3(BJ) energies (in eV) for the initial PBE0-D3(BJ) database in (a).

Figure 3

Parity plots and error distributions comparing predicted forces (in meV/Å) for (a, c) REF-HLNNP and (b, d) TL-HLNNP with computed reference values (PBE0-D3BJ/MOLOPT-TZVP) for a test set consisting of the 200 most diverse structures along the reaction pathway. Root-mean-square error (RMSE) below 150 meV/Å (force) is considered to be good performance.

Figure 4

Free energy surfaces obtained with the GLNNP (emulating PBE-D3(BJ), top) and TL-HLNNPs (emulating PBE0-D3(BJ), bottom) for the MeOH (left) and TFE (right) systems. ΔG values (in kcal/mol) presented are relative to the basin of structure 1 (0 kcal/mol). Averages and standard deviations were computed based on four independent converged dynamic simulations.

Figure 5

2D tSNE representation of the structures captured during WT-MTD simulations in methanol. Color changes from blue to red through basins A–E represent reaction progress from 1 to Int1. Enlargements of the representative structures can be found in the SI (Figures S11–S15) or visualized interactively at https://archive.materialscloud.org/record/2024.135.

Figure 6

2D tSNE representation of structures captured during WT-MTD simulations in 2,2,2-trifluoroethanol (TFE). Color changes from blue to red through basins A–E represent reaction progress from 1 to Int1. Enlargements of the representative structures can be found in the SI (Figures S11–S15) or visualized interactively at https://archive.materialscloud.org/record/2024.135.

Figure 7

Probability distribution of the number of solvent molecules explicitly involved in the proton transfer network with phthalimide (blue is MeOH and green is TFE). Statistics taken from an ensemble of 100 independent 6 ns WT-MTD, resulting in 200 TS structures for each solvent. Structures close to the TS region in TFE are provided for the two most probable TSs (with 3 and 4 solvent molecules) on the right.